Chapter 4 Diversity analysis

load("data/08032025data.Rdata")
load("data/04032025gifts.Rdata")
load("data/beta_08032025.Rdata")
genome_counts_filt <- genome_counts_filt %>% 
  column_to_rownames(var = "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0)%>% 
  select(where(~ sum(.) > 0)) %>% 
  select(-AC85) %>% 
  rownames_to_column(., "genome")

genome_tree <- keep.tip(genome_tree, tip=genome_counts_filt$genome)

table <- genome_counts_filt %>%
  t() %>%
  row_to_names(row_number = 1) %>% 
  as.data.frame() %>% 
  mutate_if(is.character,as.numeric) %>% 
  rownames_to_column(., "sample")
  
master_table <- sample_metadata %>%
  mutate(sample=Tube_code) %>%
  mutate(Tube_code=str_remove_all(Tube_code, "_a")) %>% 
  filter(Tube_code %in% table$sample) %>%
  mutate(treatment = sub("^\\d+_", "", time_point))%>% 
  left_join(., table, by=join_by("Tube_code" =="sample"))

sample_metadata <- master_table %>% 
  select(1:13)

genome_counts_filt <- master_table %>% 
  select(12,14:323) %>% 
  t() %>%
  row_to_names(row_number = 1) %>% 
  as.data.frame() %>% 
  mutate_if(is.character,as.numeric) %>% 
  filter(rowSums(. != 0, na.rm = TRUE) > 0)%>%
  dplyr::select(where(~ !all(. == 0)))%>% 
  rownames_to_column(., "genome")

genome_counts_filtering <- master_table %>% 
  select(12,14:323) %>% 
  t() %>%
  row_to_names(row_number = 1) %>% 
  as.data.frame() %>% 
  mutate_if(is.character,as.numeric) %>% 
  filter(rowSums(. != 0, na.rm = TRUE) > 0)%>%
  dplyr::select(where(~ !all(. == 0)))

genome_tree <- keep.tip(genome_tree, tip=genome_counts_filt$genome)
#match_data(genome_counts_filtering,genome_tree)

genome_metadata <- genome_metadata %>% 
  filter(genome %in% genome_counts_filt$genome)

4.1 Calculate diversities

4.1.1 Alpha diversity

# Calculate Hill numbers
richness <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 0) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(richness = 1) %>%
  rownames_to_column(var = "sample")

neutral <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(neutral = 1) %>%
  rownames_to_column(var = "sample")

phylogenetic <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, tree = genome_tree) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(phylogenetic = 1) %>%
  rownames_to_column(var = "sample")

genome_gifts1 <- genome_gifts[genome_counts_filt$genome,]
genome_gifts1 <- genome_gifts1[, colSums(genome_gifts1 != 0) > 0]

genome_counts_filt1 <- genome_counts_filt[genome_counts_filt$genome %in% rownames(genome_gifts1),]
genome_counts_filt1 <- genome_counts_filt1 %>% 
  remove_rownames() %>% 
  column_to_rownames(var = "genome") %>% 
      filter(rowSums(. != 0, na.rm = TRUE) > 0) %>% 
  rownames_to_column(., "genome")

# Merge all metrics
alpha_div <- richness %>%
  full_join(neutral, by = join_by(sample == sample)) %>%
  full_join(phylogenetic, by = join_by(sample == sample))
save(alpha_div, 
     file = "data/alpha_08032025.Rdata")

4.1.2 Beta diversity

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, tree = genome_tree)

genome_gifts1 <- genome_gifts[genome_counts_filt$genome,]
genome_gifts1 <- genome_gifts1[, colSums(genome_gifts1 != 0) > 0]
genome_counts_filt1 <- genome_counts_filt[genome_counts_filt$genome %in% rownames(genome_gifts1),]
save(beta_q0n, 
     beta_q1n, 
     beta_q1p, 
     file = "data/beta_08032025.Rdata")

4.2 Does captivity change the microbial community?

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  filter(time_point %in% c("Acclimation", "Wild") ) 

4.2.1 Alpha diversity

label_map <- c(
  "Cold_wet" = "Cold population", 
  "Hot_dry" = "Warm-population",
  "richness" = "Species Richness",
  "neutral" = "Neutral Diversity",
  "phylogenetic" = "Phylogenetic Diversity")

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(
  metric = factor(metric, levels = c("richness", "neutral", "phylogenetic")),
  time_point = factor(time_point, levels = c("Wild", "Acclimation"))) %>%
  filter(time_point %in% c("Wild", "Acclimation")) %>% 
  ggplot(aes(y = value, x = time_point, color=Population, fill=Population)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(alpha=0.5) +
  scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
  facet_grid(metric ~ Population, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
  theme_classic() +
  theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

Richness

Modelq0GLMMNB <- glmer.nb(richness ~ Population*time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(9.4635)  ( log )
Formula: richness ~ Population * time_point + (1 | individual)
   Data: alpha_div_meta

     AIC      BIC   logLik deviance df.resid 
   518.4    530.2   -253.2    506.4       47 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.07994 -0.34567  0.00022  0.55317  1.22959 

Random effects:
 Groups     Name        Variance Std.Dev.
 individual (Intercept) 0.1126   0.3355  
Number of obs: 53, groups:  individual, 27

Fixed effects:
                                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)                        3.7707     0.1208  31.218  < 2e-16 ***
PopulationHot_dry                  0.6756     0.2006   3.368 0.000757 ***
time_pointWild                     0.4312     0.1223   3.525 0.000424 ***
PopulationHot_dry:time_pointWild  -0.5331     0.2048  -2.604 0.009221 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) PpltH_ tm_pnW
PpltnHt_dry -0.602              
time_pntWld -0.557  0.335       
PpltnHt_:_W  0.338 -0.526 -0.601
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ time_point)
$emmeans
 time_point  emmean    SE  df asymp.LCL asymp.UCL
 Acclimation   4.11 0.100 Inf      3.91      4.31
 Wild          4.27 0.099 Inf      4.08      4.47

Results are averaged over the levels of: Population 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast           estimate    SE  df z.ratio p.value
 Acclimation - Wild   -0.165 0.102 Inf  -1.615  0.1063

Results are averaged over the levels of: Population 
Results are given on the log (not the response) scale. 
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ Population)
$emmeans
 Population emmean    SE  df asymp.LCL asymp.UCL
 Cold_wet     3.99 0.101 Inf      3.79      4.18
 Hot_dry      4.40 0.138 Inf      4.12      4.67

Results are averaged over the levels of: time_point 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast           estimate    SE  df z.ratio p.value
 Cold_wet - Hot_dry   -0.409 0.171 Inf  -2.396  0.0166

Results are averaged over the levels of: time_point 
Results are given on the log (not the response) scale. 

Neutral

Model_neutral <- nlme::lme(fixed = neutral ~ time_point*Population, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_neutral)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  402.7832 414.1342 -195.3916

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:     5.51778 10.54924

Fixed effects:  neutral ~ time_point * Population 
                                     Value Std.Error DF   t-value p-value
(Intercept)                       19.37729  2.883717 25  6.719554  0.0000
time_pointWild                    17.27881  3.578682 24  4.828260  0.0001
PopulationHot_dry                 24.71329  4.905493 25  5.037882  0.0000
time_pointWild:PopulationHot_dry -22.28917  6.126768 24 -3.637998  0.0013
 Correlation: 
                                 (Intr) tm_pnW PpltH_
time_pointWild                   -0.642              
PopulationHot_dry                -0.588  0.377       
time_pointWild:PopulationHot_dry  0.375 -0.584 -0.632

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-2.5974563 -0.5606149  0.1376820  0.5636829  1.6703831 

Number of Observations: 53
Number of Groups: 27 
emmeans::emmeans(Model_neutral, pairwise ~ time_point)
$emmeans
 time_point  emmean   SE df lower.CL upper.CL
 Acclimation   31.7 2.45 25     26.7     36.8
 Wild          37.9 2.43 24     32.9     42.9

Results are averaged over the levels of: Population 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast           estimate   SE df t.ratio p.value
 Acclimation - Wild    -6.13 3.06 24  -2.002  0.0567

Results are averaged over the levels of: Population 
Degrees-of-freedom method: containment 
emmeans::emmeans(Model_neutral, pairwise ~ Population)
$emmeans
 Population emmean   SE df lower.CL upper.CL
 Cold_wet     28.0 2.21 24     23.5     32.6
 Hot_dry      41.6 3.09 24     35.2     48.0

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast           estimate  SE df t.ratio p.value
 Cold_wet - Hot_dry    -13.6 3.8 24  -3.568  0.0016

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 

Phylogenetic

Model_phylo <- nlme::lme(fixed = phylogenetic ~ time_point*Population, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_phylo)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  193.7567 205.1077 -90.87837

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:   0.7860836 1.191838

Fixed effects:  phylogenetic ~ time_point * Population 
                                     Value Std.Error DF   t-value p-value
(Intercept)                       5.343318 0.3453897 25 15.470404  0.0000
time_pointWild                   -0.135493 0.4048212 24 -0.334698  0.7408
PopulationHot_dry                 1.119702 0.5880336 25  1.904146  0.0685
time_pointWild:PopulationHot_dry -0.602016 0.6924897 24 -0.869350  0.3933
 Correlation: 
                                 (Intr) tm_pnW PpltH_
time_pointWild                   -0.608              
PopulationHot_dry                -0.587  0.357       
time_pointWild:PopulationHot_dry  0.355 -0.585 -0.596

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-2.2027374 -0.4577650  0.0610368  0.3504732  1.7675077 

Number of Observations: 53
Number of Groups: 27 
emmeans::emmeans(Model_neutral, pairwise ~ time_point)
$emmeans
 time_point  emmean   SE df lower.CL upper.CL
 Acclimation   31.7 2.45 25     26.7     36.8
 Wild          37.9 2.43 24     32.9     42.9

Results are averaged over the levels of: Population 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast           estimate   SE df t.ratio p.value
 Acclimation - Wild    -6.13 3.06 24  -2.002  0.0567

Results are averaged over the levels of: Population 
Degrees-of-freedom method: containment 
emmeans::emmeans(Model_neutral, pairwise ~ Population)
$emmeans
 Population emmean   SE df lower.CL upper.CL
 Cold_wet     28.0 2.21 24     23.5     32.6
 Hot_dry      41.6 3.09 24     35.2     48.0

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast           estimate  SE df t.ratio p.value
 Cold_wet - Hot_dry    -13.6 3.8 24  -3.568  0.0016

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 

4.2.2 Beta diversity

samples_to_keep <- sample_metadata %>%
  filter(time_point %in% c("Acclimation", "Wild")) %>% 
  select(Tube_code) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(time_point %in% c("Acclimation", "Wild"))

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.03445 0.034446 4.6041    999  0.036 *
Residuals 51 0.38156 0.007482                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(richness ~ time_point*Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
time_point 1 0.5784094 0.03513500 2.236493 0.001
Population 1 2.8737192 0.17456169 11.111600 0.001
time_point:Population 1 0.3378127 0.02052015 1.306196 0.001
Residual 49 12.6725436 0.76978316 NA NA
Total 52 16.4624849 1.00000000 NA NA
subset_meta_arrange <- column_to_rownames(subset_meta, "Tube_code")
subset_meta_arrange<-subset_meta_arrange[labels(richness),]

pairwise <- pairwise.adonis(richness, subset_meta_arrange$Population, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Cold_wet vs Hot_dry 1 2.87727 10.8015 0.1747774 0.001 0.001 **
pairwise <- pairwise.adonis(richness, subset_meta_arrange$time_point, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Wild vs Acclimation 1 0.5784094 1.857135 0.035135 0.023 0.023 .

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$time_point) %>% permutest(.)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01130 0.0113018 1.2646    999  0.259
Residuals 51 0.45578 0.0089369                     
adonis2(neutral ~ time_point*Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
time_point 1 0.8535852 0.05306573 3.688514 0.001
Population 1 3.3711165 0.20957575 14.567274 0.001
time_point:Population 1 0.5212922 0.03240772 2.252609 0.001
Residual 49 11.3394385 0.70495080 NA NA
Total 52 16.0854325 1.00000000 NA NA
subset_meta_arrange <- column_to_rownames(subset_meta, "Tube_code")
subset_meta_arrange<-subset_meta_arrange[labels(neutral),]

pairwise <- pairwise.adonis(neutral, subset_meta_arrange$Population, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Cold_wet vs Hot_dry 1 3.369473 13.51397 0.2094736 0.001 0.001 **
pairwise <- pairwise.adonis(neutral, subset_meta_arrange$time_point, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Wild vs Acclimation 1 0.8535852 2.858015 0.05306573 0.002 0.002 *

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq     F N.Perm Pr(>F)
Groups     1 0.00207 0.0020675 0.125    999  0.735
Residuals 51 0.84374 0.0165440                    
adonis2(phylo ~ time_point*Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
time_point 1 0.2282721 0.07008479 4.486412 0.001
Population 1 0.3478468 0.10679697 6.836508 0.001
time_point:Population 1 0.1878081 0.05766139 3.691140 0.003
Residual 49 2.4931580 0.76545685 NA NA
Total 52 3.2570850 1.00000000 NA NA
subset_meta_arrange <- column_to_rownames(subset_meta, "Tube_code")
subset_meta_arrange<-subset_meta_arrange[labels(phylo),]

pairwise <- pairwise.adonis(phylo, subset_meta_arrange$Population, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Cold_wet vs Hot_dry 1 0.3462509 6.066576 0.106307 0.002 0.002 *
pairwise <- pairwise.adonis(phylo, subset_meta_arrange$time_point, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Wild vs Acclimation 1 0.2282721 3.84371 0.07008479 0.004 0.004 *

dbRDA

#Richness
cca_ord <- capscale(formula = richness ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                 "subset_meta$time_pointWild" = "Wild",
                                 "subset_meta$PopulationHot_dry"="Warm population",
                                 "subset_meta$time_pointWild:subset_meta$PopulationHot_dry"="Interaction Wam population")

beta_richness_rda <-CAP_df %>%
  group_by(Population, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
  scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()

#Neutral
cca_ord <- capscale(formula = neutral ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                 "subset_meta$time_pointWild" = "Wild",
                                 "subset_meta$PopulationHot_dry"="Warm population",
                                 "subset_meta$time_pointWild:subset_meta$PopulationHot_dry"="Interaction Wam population")

beta_neutral_rda <- CAP_df %>%
  group_by(Population, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
  scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50"))+
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()


#Phylogenetic
cca_ord <- capscale(formula = phylo ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                 "subset_meta$time_pointWild" = "Wild",
                                 "subset_meta$PopulationHot_dry"="Warm population",
                                 "subset_meta$time_pointWild:subset_meta$PopulationHot_dry"="Interaction Wam population")

beta_phylogenetic_rda<- CAP_df %>%
  group_by(type, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
  scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()
ggarrange(beta_richness_rda, beta_neutral_rda, beta_phylogenetic_rda, ncol=3, nrow=1, common.legend = TRUE, legend="right")

4.2.3 Cold population

4.2.3.1 Structural zeros

wild_samples <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point == "Wild") %>% 
  dplyr::select(sample) %>%
  pull()

acclimation_samples <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point == "Acclimation") %>% 
  dplyr::select(sample) %>% pull()

structural_zeros <- genome_counts_filt %>% 
  select(one_of(c("genome",subset_meta$sample))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>% 
   rowwise() %>% #compute for each row (genome)
   mutate(all_zeros_wild = all(c_across(all_of(wild_samples)) == 0)) %>% 
   mutate(all_zeros_acclimation= all(c_across(all_of(acclimation_samples)) == 0)) %>% 
   mutate(average_wild = mean(c_across(all_of(wild_samples)), na.rm = TRUE)) %>% 
   mutate(average_acclimation = mean(c_across(all_of(acclimation_samples)), na.rm = TRUE)) %>% 
   filter(all_zeros_wild == TRUE || all_zeros_acclimation==TRUE)  %>% 
   mutate(present = case_when(
      all_zeros_wild & !all_zeros_acclimation ~ "acclimation",
      !all_zeros_wild & all_zeros_acclimation ~ "wild",
      !all_zeros_wild & !all_zeros_acclimation ~ "None",
      TRUE ~ NA_character_
    )) %>%
   mutate(average = ifelse(present == "acclimation", average_wild, average_acclimation)) %>%
   dplyr::select(genome, present, average) %>%
   left_join(genome_metadata, by=join_by(genome==genome)) %>%
   arrange(present,-average)

Wild

structural_zeros %>% 
  filter(present=="wild") %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count)) 
# A tibble: 5 × 2
# Rowwise: 
  phylum           sample_count
  <chr>                   <int>
1 p__Bacillota_A              7
2 p__Bacillota                1
3 p__Bacillota_B              1
4 p__Bacteroidota             1
5 p__Spirochaetota            1

Acclimation

structural_zeros %>% 
  filter(present=="acclimation") %>% 
  select(1,5:10)
# A tibble: 14 × 7
# Rowwise: 
   genome                phylum            class                  order               family                genus              species                       
   <chr>                 <chr>             <chr>                  <chr>               <chr>                 <chr>              <chr>                         
 1 AH1_2nd_12:bin_000001 p__Pseudomonadota c__Gammaproteobacteria o__Pseudomonadales  f__Moraxellaceae      g__Acinetobacter   s__Acinetobacter johnsonii    
 2 AH1_2nd_15:bin_000001 p__Pseudomonadota c__Alphaproteobacteria o__Rhizobiales      f__Rhizobiaceae       g__Agrobacterium   s__Agrobacterium tumefaciens_H
 3 AH1_2nd_17:bin_000010 p__Bacteroidota   c__Bacteroidia         o__Bacteroidales    f__Rikenellaceae      g__Mucinivorans    s__                           
 4 AH1_2nd_17:bin_000038 p__Bacteroidota   c__Bacteroidia         o__Bacteroidales    f__Bacteroidaceae     g__Bacteroides_G   s__                           
 5 AH1_2nd_20:bin_000014 p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales f__Enterobacteriaceae g__Citrobacter     s__Citrobacter portucalensis  
 6 AH1_2nd_20:bin_000061 p__Bacillota      c__Bacilli             o__Lactobacillales  f__Enterococcaceae    g__Enterococcus    s__                           
 7 AH1_2nd_20:bin_000073 p__Bacillota      c__Bacilli             o__Lactobacillales  f__Enterococcaceae    g__Enterococcus    s__Enterococcus sp002174455   
 8 LI1_2nd_10:bin_000001 p__Bacillota      c__Bacilli             o__Lactobacillales  f__Enterococcaceae    g__Enterococcus_A  s__Enterococcus_A raffinosus  
 9 LI1_2nd_10:bin_000017 p__Chlamydiota    c__Chlamydiia          o__Chlamydiales     f__Chlamydiaceae      g__                s__                           
10 LI1_2nd_2:bin_000039  p__Bacillota_A    c__Clostridia          o__TANB77           f__CAG-508            g__                s__                           
11 LI1_2nd_7:bin_000045  p__Bacillota_A    c__Clostridia          o__Oscillospirales  f__Acutalibacteraceae g__                s__                           
12 LI1_2nd_7:bin_000074  p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales f__Enterobacteriaceae g__Escherichia     s__Escherichia coli           
13 LI1_2nd_8:bin_000030  p__Actinomycetota c__Actinomycetia       o__Mycobacteriales  f__Mycobacteriaceae   g__Corynebacterium s__                           
14 LI1_2nd_8:bin_000079  p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales f__Enterobacteriaceae g__Citrobacter_A   s__Citrobacter_A amalonaticus 

Community elements differences in structural zeros

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))
uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)
# A tibble: 0 × 3
# ℹ 3 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>

4.2.3.2 Differential abundance analysis: composition

phylo_samples <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point %in% c("Wild", "Acclimation") )%>% 
  column_to_rownames("sample") %>% 
  sample_data()
phylo_genome <- genome_counts_filt %>% 
  filter(!genome %in% structural_zeros$genome) %>% 
  select(one_of(c("genome",rownames(phylo_samples)))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
  column_to_rownames("genome") %>% 
  otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>% 
  filter(genome %in% rownames(phylo_genome)) %>% 
  column_to_rownames("genome") %>% 
  dplyr::select(domain,phylum,class,order,family,genus,species) %>% 
  as.matrix() %>% 
  tax_table()

physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)
ancom_rand_output = ancombc2(data = physeq_genome_filtered, 
                  assay_name = "counts", 
                  tax_level = NULL,
                  fix_formula = "time_point",
                  rand_formula = "(1|individual)",
                  p_adj_method = "holm", 
                  pseudo_sens = TRUE,
                  prv_cut =0, 
                  lib_cut = 0, 
                  s0_perc = 0.05,
                  group = NULL, 
                  struc_zero = FALSE, 
                  neg_lb = FALSE,
                  alpha = 0.05, 
                  n_cl = 2, 
                  verbose = TRUE,
                  global = FALSE, 
                  pairwise = FALSE, 
                  dunnet = FALSE, 
                  trend = FALSE,
                  iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
                  em_control = list(tol = 1e-5, max_iter = 100),
                  lme_control = lme4::lmerControl(),
                  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100), 
                  trend_control = NULL)
save(ancom_rand_output, 
     file = "data/ancombc_wild_accl.Rdata")
load("data/ancombc_wild_accl.Rdata")
taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
  rownames_to_column(., "taxon") %>%
  mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))

ancombc_rand_table_mag <- ancom_rand_output$res %>%
  dplyr::select(taxon, lfc_time_point1_Acclimation, p_time_point1_Acclimation) %>%
  filter(p_time_point1_Acclimation < 0.05) %>%
  dplyr::arrange(p_time_point1_Acclimation) %>%
  merge(., taxonomy, by="taxon") %>%
  mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
  dplyr::arrange(lfc_time_point1_Acclimation)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", ""))  %>%
  right_join(taxonomy, by=join_by(phylum == phylum)) %>%
  dplyr::select(phylum, colors) %>%
  mutate(colors = str_c(colors, "80"))  %>% #add 80% alpha
    unique() %>%
    dplyr::arrange(phylum)

tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
  
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
    dplyr::arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()
ancombc_rand_table_mag%>%
      mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_point1_Acclimation, y=forcats::fct_reorder(genome,lfc_time_point1_Acclimation), fill=phylum)) + #forcats::fct_rev()
  geom_col() + 
  scale_fill_manual(values=tax_color) + 
  geom_hline(yintercept=0) + 
#  coord_flip()+
  theme(panel.background = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 8),
        axis.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14, face = "bold"),
        legend.position = "right", legend.box = "vertical")+
  xlab("log2FoldChange") + 
  ylab("Species")+
  guides(fill=guide_legend(title="Phylum"))

Phyla of the significant MAGs

ancombc_rand_table_mag%>%
  filter(lfc_time_point1_Acclimation<0)  %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count))  
        phylum sample_count
1  Bacillota_A           16
2 Bacteroidota            9
3  Bacillota_C            1
ancombc_rand_table_mag%>%
  filter(lfc_time_point1_Acclimation<0)  %>% 
  select(phylum, genus, species)
         phylum                genus                      species
1   Bacillota_A           MGBC140009                             
2   Bacillota_A                                                  
3   Bacillota_A                                                  
4   Bacillota_A     Negativibacillus                             
5  Bacteroidota      Parabacteroides                             
6  Bacteroidota          Bacteroides                             
7   Bacillota_A           Hungatella                             
8  Bacteroidota          Bacteroides                             
9  Bacteroidota      Parabacteroides                             
10 Bacteroidota      Parabacteroides                             
11 Bacteroidota          Bacteroides                             
12  Bacillota_A                                                  
13 Bacteroidota          Bacteroides                             
14  Bacillota_A                                                  
15  Bacillota_A           Copromonas                             
16  Bacillota_A        Enterocloster                             
17  Bacillota_A      Velocimicrobium                             
18  Bacillota_A               CAG-95                             
19  Bacillota_C                                                  
20 Bacteroidota          Bacteroides                             
21  Bacillota_A                                                  
22  Bacillota_A           Copromonas                             
23 Bacteroidota          Bacteroides Bacteroides thetaiotaomicron
24  Bacillota_A Pseudoflavonifractor                             
25  Bacillota_A        Enterocloster                             
26  Bacillota_A             JALFVM01                             
ancombc_rand_table_mag%>%
  filter(lfc_time_point1_Acclimation<0)  %>% 
  count(genus, name = "sample_count") %>%
  arrange(desc(sample_count))
                  genus sample_count
1                                  6
2           Bacteroides            6
3       Parabacteroides            3
4            Copromonas            2
5         Enterocloster            2
6                CAG-95            1
7            Hungatella            1
8              JALFVM01            1
9            MGBC140009            1
10     Negativibacillus            1
11 Pseudoflavonifractor            1
12      Velocimicrobium            1

4.2.3.3 Differences in the functional capacity

GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>% 
  select_if(~ !is.numeric(.) || sum(.) != 0)

GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)

GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)

sample_sub <- sample_metadata %>%
  filter(Population == "Cold_wet" & treatment %in% c("Acclimation", "Wild"))

genome_counts_row <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% 
    select(one_of(c("genome",sample_sub$sample))) %>%
  column_to_rownames(., "genome") 

GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)
4.2.3.3.1 Community elements differences:

MCI

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))
# A tibble: 2 × 3
  treatment     MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.335 0.0324
2 Wild        0.350 0.0237
MCI <- GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.94347, p-value = 0.07144
wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 114, p-value = 0.2066
alternative hypothesis: true location shift is not equal to 0
element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))
uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

significant_elements <- element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)%>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "sample")%>% 
  left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))

difference_table <- element_gift_filt %>%
  select(-sample) %>%
  group_by(treatment) %>%
  summarise(across(everything(), mean)) %>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric) %>%
  rownames_to_column(., "Elements") %>%
  left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>% 
  arrange(Function) %>% 
  mutate(Difference=Acclimation-Wild)%>% 
  mutate(group_color = ifelse(Difference <0, "Wild","Acclimation")) 
difference_table %>%
  ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) + 
  geom_col() +
#  geom_point(size=4) + 
  scale_fill_manual(values=c('#008080', '#00808050')) +
  geom_hline(yintercept=0) + 
  coord_flip()+
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        legend.position = "right", 
        panel.background = element_blank(),
          panel.grid.major = element_line(size = 0.15, linetype = 'solid',
                                colour = "grey"))+
  xlab("Microbial Functional Capacity") + 
  ylab("Mean difference")+
  labs(fill="Treatment")

4.2.3.3.2 Community functions differences

MCI

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))
# A tibble: 2 × 3
  treatment     MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.330 0.0320
2 Wild        0.346 0.0196
MCI <- GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.94507, p-value = 0.07998
wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 100, p-value = 0.08313
alternative hypothesis: true location shift is not equal to 0
function_gift <- GIFTs_functions_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  merge(., sample_metadata[c(12,13)], by="sample")
unique_funct_db<- GIFT_db[c(3,4,5)] %>% 
  distinct(Code_function, .keep_all = TRUE)

significant_functional <- function_gift %>%
  pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
  mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
  filter(p_adjust < 0.05)%>%
  left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))
Code_function Acclimation Wild Function Difference treatment
B03 0.1097356 0.1256091 Amino acid derivative biosynthesis -0.0158735 Wild
D03 0.2921075 0.3442827 Sugar degradation -0.0521752 Wild

4.2.4 Warm population

4.2.4.1 Structural zeros

wild_samples_warm <- sample_metadata %>%
  filter(Population == "Hot_dry" & time_point == "Wild") %>% 
  dplyr::select(sample) %>%
  pull()

acclimation_samples_warm <- sample_metadata %>%
  filter(Population == "Hot_dry" & time_point == "Acclimation") %>% 
  dplyr::select(sample) %>% pull()

structural_zeros <- genome_counts_filt %>% 
  select(one_of(c("genome",subset_meta$sample))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>% 
   rowwise() %>% #compute for each row (genome)
   mutate(all_zeros_wild_warm = all(c_across(all_of(wild_samples_warm)) == 0)) %>% 
   mutate(all_zeros_acclimation_warm= all(c_across(all_of(acclimation_samples_warm)) == 0)) %>% 
   mutate(average_wild_warm = mean(c_across(all_of(wild_samples_warm)), na.rm = TRUE)) %>% 
   mutate(average_acclimation_warm = mean(c_across(all_of(acclimation_samples_warm)), na.rm = TRUE)) %>% 
   filter(all_zeros_wild_warm == TRUE || all_zeros_acclimation_warm==TRUE)  %>% 
   mutate(present = case_when(
      all_zeros_wild_warm & !all_zeros_acclimation_warm ~ "acclimation",
      !all_zeros_wild_warm & all_zeros_acclimation_warm ~ "wild",
      !all_zeros_wild_warm & !all_zeros_acclimation_warm ~ "None",
      TRUE ~ NA_character_
    )) %>%
   mutate(average = ifelse(present == "acclimation", average_wild_warm, average_acclimation_warm)) %>%
   dplyr::select(genome, present, average) %>%
   left_join(genome_metadata, by=join_by(genome==genome)) %>%
   arrange(present,-average)

Wild

structural_zeros %>% 
  filter(present=="wild") %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count)) 
# A tibble: 7 × 2
# Rowwise: 
  phylum              sample_count
  <chr>                      <int>
1 p__Bacillota_A                11
2 p__Bacillota                   4
3 p__Pseudomonadota              4
4 p__Bacillota_C                 2
5 p__Bacteroidota                1
6 p__Desulfobacterota            1
7 p__Fusobacteriota              1

Acclimation

structural_zeros %>% 
  filter(present=="acclimation") %>% 
  select(1,5:10)
# A tibble: 32 × 7
# Rowwise: 
   genome                phylum             class                  order                  family                  genus               species                       
   <chr>                 <chr>              <chr>                  <chr>                  <chr>                   <chr>               <chr>                         
 1 AH1_2nd_10:bin_000010 p__Cyanobacteriota c__Vampirovibrionia    o__Gastranaerophilales f__Gastranaerophilaceae g__CALURL01         s__                           
 2 AH1_2nd_12:bin_000001 p__Pseudomonadota  c__Gammaproteobacteria o__Pseudomonadales     f__Moraxellaceae        g__Acinetobacter    s__Acinetobacter johnsonii    
 3 AH1_2nd_14:bin_000015 p__Bacillota_A     c__Clostridia          o__Christensenellales  f__UBA3700              g__                 s__                           
 4 AH1_2nd_14:bin_000063 p__Bacillota_A     c__Clostridia          o__Lachnospirales      f__Lachnospiraceae      g__Blautia_A        s__                           
 5 AH1_2nd_15:bin_000001 p__Pseudomonadota  c__Alphaproteobacteria o__Rhizobiales         f__Rhizobiaceae         g__Agrobacterium    s__Agrobacterium tumefaciens_H
 6 AH1_2nd_16:bin_000033 p__Bacteroidota    c__Bacteroidia         o__Bacteroidales       f__                     g__                 s__                           
 7 AH1_2nd_16:bin_000096 p__Bacillota_A     c__Clostridia          o__Oscillospirales     f__Ruminococcaceae      g__Ruthenibacterium s__                           
 8 AH1_2nd_17:bin_000015 p__Bacillota       c__Bacilli             o__Erysipelotrichales  f__Erysipelotrichaceae  g__                 s__                           
 9 AH1_2nd_17:bin_000023 p__Bacillota       c__Bacilli             o__Erysipelotrichales  f__Coprobacillaceae     g__Coprobacillus    s__                           
10 AH1_2nd_18:bin_000039 p__Bacteroidota    c__Bacteroidia         o__Bacteroidales       f__Rikenellaceae        g__Alistipes        s__                           
# ℹ 22 more rows

Community elements differences in structural zeros

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))
uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)
# A tibble: 0 × 3
# ℹ 3 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>

4.2.5 Differential abundance analysis: composition

phylo_samples_warm <- sample_metadata %>%
  filter(Population == "Hot_dry" & time_point %in% c("Wild", "Acclimation") )%>% 
  column_to_rownames("sample") %>% 
  sample_data()
phylo_genome_warm <- genome_counts_filt %>% 
  filter(!genome %in% structural_zeros$genome) %>% 
  select(one_of(c("genome",rownames(phylo_samples_warm)))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
  column_to_rownames("genome") %>% 
  otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy_warm <- genome_metadata %>% 
  filter(genome %in% rownames(phylo_genome_warm)) %>% 
  column_to_rownames("genome") %>% 
  dplyr::select(domain,phylum,class,order,family,genus,species) %>% 
  as.matrix() %>% 
  tax_table()

physeq_genome_filtered_warm <- phyloseq(phylo_genome_warm, phylo_taxonomy_warm, phylo_samples_warm)
ancom_rand_output_warm = ancombc2(data = physeq_genome_filtered_warm, 
                  assay_name = "counts", 
                  tax_level = NULL,
                  fix_formula = "time_point",
                  rand_formula = "(1|individual)",
                  p_adj_method = "holm", 
                  pseudo_sens = TRUE,
                  prv_cut =0, 
                  lib_cut = 0, 
                  s0_perc = 0.05,
                  group = NULL, 
                  struc_zero = FALSE, 
                  neg_lb = FALSE,
                  alpha = 0.05, 
                  n_cl = 2, 
                  verbose = TRUE,
                  global = FALSE, 
                  pairwise = FALSE, 
                  dunnet = FALSE, 
                  trend = FALSE,
                  iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
                  em_control = list(tol = 1e-5, max_iter = 100),
                  lme_control = lme4::lmerControl(),
                  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100), 
                  trend_control = NULL)
save(ancom_rand_output_warm, 
     file = "data/ancombc_wild_accl_warm.Rdata")
load("data/ancombc_wild_accl_warm.Rdata")
taxonomy <- data.frame(physeq_genome_filtered_warm@tax_table) %>%
  rownames_to_column(., "taxon") %>%
  mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))

ancombc_rand_table_mag <- ancom_rand_output_warm$res %>%
  dplyr::select(taxon, lfc_time_point1_Acclimation, p_time_point1_Acclimation) %>%
  filter(p_time_point1_Acclimation < 0.05) %>%
  dplyr::arrange(p_time_point1_Acclimation) %>%
  merge(., taxonomy, by="taxon") %>%
  mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
  dplyr::arrange(lfc_time_point1_Acclimation)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", ""))  %>%
  right_join(taxonomy, by=join_by(phylum == phylum)) %>%
  dplyr::select(phylum, colors) %>%
  mutate(colors = str_c(colors, "80"))  %>% #add 80% alpha
    unique() %>%
    dplyr::arrange(phylum)

tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
  
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
    dplyr::arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()
ancombc_rand_table_mag%>%
      mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_point1_Acclimation, y=forcats::fct_reorder(genome,lfc_time_point1_Acclimation), fill=phylum)) + #forcats::fct_rev()
  geom_col() + 
  scale_fill_manual(values=tax_color) + 
  geom_hline(yintercept=0) + 
#  coord_flip()+
  theme(panel.background = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 8),
        axis.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14, face = "bold"),
        legend.position = "right", legend.box = "vertical")+
  xlab("log2FoldChange") + 
  ylab("Species")+
  guides(fill=guide_legend(title="Phylum"))

Phyla of the significant MAGs

ancombc_rand_table_mag%>%
  filter(lfc_time_point1_Acclimation<0)  %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count))  

ancombc_rand_table_mag%>%
  filter(lfc_time_point1_Acclimation<0)  %>% 
  select(phylum, genus, species)

ancombc_rand_table_mag%>%
  filter(lfc_time_point1_Acclimation<0)  %>% 
  count(genus, name = "sample_count") %>%
  arrange(desc(sample_count))

4.2.6 Differences in the functional capacity

GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>% 
  select_if(~ !is.numeric(.) || sum(.) != 0)

GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)

GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)

sample_sub <- sample_metadata %>%
  filter(Population == "Hot_dry" & treatment %in% c("Acclimation", "Wild"))

genome_counts_row <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% 
    select(one_of(c("genome",sample_sub$sample))) %>%
  column_to_rownames(., "genome") 

GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)

4.2.6.1 Community elements differences:

MCI

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))
# A tibble: 2 × 3
  treatment     MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.355 0.0228
2 Wild        0.317 0.0335
MCI <- GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.94062, p-value = 0.297
wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 68, p-value = 0.01419
alternative hypothesis: true location shift is not equal to 0
element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))
uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

significant_elements <- element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)%>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "sample")%>% 
  left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))

difference_table <- element_gift_filt %>%
  select(-sample) %>%
  group_by(treatment) %>%
  summarise(across(everything(), mean)) %>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric) %>%
  rownames_to_column(., "Elements") %>%
  left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>% 
  arrange(Function) %>% 
  mutate(Difference=Acclimation-Wild)%>% 
  mutate(group_color = ifelse(Difference <0, "Wild","Acclimation")) 
difference_table %>%
  ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) + 
  geom_col() +
#  geom_point(size=4) + 
  scale_fill_manual(values=c("#d57d2c50","#d57d2c")) +
  geom_hline(yintercept=0) + 
  coord_flip()+
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        legend.position = "right", 
        panel.background = element_blank(),
          panel.grid.major = element_line(size = 0.15, linetype = 'solid',
                                colour = "grey"))+
  xlab("Microbial Functional Capacity") + 
  ylab("Mean difference")+
  labs(fill="Treatment")

4.2.6.2 Community functions differences

MCI

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))
# A tibble: 2 × 3
  treatment     MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.348 0.0158
2 Wild        0.327 0.0245
MCI <- GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.95485, p-value = 0.5061
wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 61, p-value = 0.07701
alternative hypothesis: true location shift is not equal to 0
function_gift <- GIFTs_functions_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  merge(., sample_metadata[c(12,13)], by="sample")
unique_funct_db<- GIFT_db[c(3,4,5)] %>% 
  distinct(Code_function, .keep_all = TRUE)

significant_functional <- function_gift %>%
  pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
  mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
  filter(p_adjust < 0.05)%>%
  left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))
Code_function Acclimation Wild Function Difference treatment
S02 0.1699218 0.2472264 Appendages -0.0773046 Wild

4.3 Does the antibiotic administration change the microbial community?

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  filter(time_point %in% c("Acclimation", "Antibiotics") ) 

4.3.1 Alpha diversity

label_map <- c(
  "Cold_wet" = "Cold population", 
  "Hot_dry" = "Warm-population",
  "richness" = "Species Richness",
  "neutral" = "Neutral Diversity",
  "phylogenetic" = "Phylogenetic Diversity")

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
  filter(time_point %in% c("Acclimation", "Antibiotics") ) %>% 
  ggplot(aes(y = value, x = time_point, color=Population, fill=Population)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(alpha=0.5) +
  scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
  facet_grid(metric ~ Population, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
  theme_classic() +
  theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

Richness

Modelq0GLMMNB <- glmer.nb(richness ~ time_point*Population+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(5.1645)  ( log )
Formula: richness ~ time_point * Population + (1 | individual)
   Data: alpha_div_meta

     AIC      BIC   logLik deviance df.resid 
   434.8    446.1   -211.4    422.8       43 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.61979 -0.56608  0.06074  0.51882  2.39923 

Random effects:
 Groups     Name        Variance Std.Dev.
 individual (Intercept) 0.2193   0.4683  
Number of obs: 49, groups:  individual, 27

Fixed effects:
                                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)                               3.7459     0.1623  23.073  < 2e-16 ***
time_pointAntibiotics                    -1.1523     0.1862  -6.190 6.02e-10 ***
PopulationHot_dry                         0.7055     0.2719   2.595  0.00946 ** 
time_pointAntibiotics:PopulationHot_dry  -0.3586     0.3037  -1.181  0.23769    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) tm_pnA PpltH_
tm_pntAntbt -0.466              
PpltnHt_dry -0.598  0.278       
tm_pntA:PH_  0.292 -0.611 -0.461
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ time_point)
$emmeans
 time_point  emmean    SE  df asymp.LCL asymp.UCL
 Acclimation   4.10 0.136 Inf      3.83      4.36
 Antibiotics   2.77 0.151 Inf      2.47      3.06

Results are averaged over the levels of: Population 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast                  estimate    SE  df z.ratio p.value
 Acclimation - Antibiotics     1.33 0.152 Inf   8.751  <.0001

Results are averaged over the levels of: Population 
Results are given on the log (not the response) scale. 
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ Population)
$emmeans
 Population emmean    SE  df asymp.LCL asymp.UCL
 Cold_wet     3.17 0.145 Inf      2.89      3.45
 Hot_dry      3.70 0.196 Inf      3.31      4.08

Results are averaged over the levels of: time_point 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast           estimate    SE  df z.ratio p.value
 Cold_wet - Hot_dry   -0.526 0.243 Inf  -2.168  0.0302

Results are averaged over the levels of: time_point 
Results are given on the log (not the response) scale. 

Neutral

Model_neutral <- nlme::lme(fixed = neutral ~ time_point*Population, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_neutral)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  342.5018 353.3418 -165.2509

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    3.282619 7.923552

Fixed effects:  neutral ~ time_point * Population 
                                            Value Std.Error DF   t-value p-value
(Intercept)                              19.11713  2.078908 25  9.195758   0e+00
time_pointAntibiotics                   -11.46317  2.832826 20 -4.046551   6e-04
PopulationHot_dry                        24.97345  3.534826 25  7.064972   0e+00
time_pointAntibiotics:PopulationHot_dry -20.91701  4.793363 20 -4.363745   3e-04
 Correlation: 
                                        (Intr) tm_pnA PpltH_
time_pointAntibiotics                   -0.633              
PopulationHot_dry                       -0.588  0.373       
time_pointAntibiotics:PopulationHot_dry  0.374 -0.591 -0.632

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.8258256 -0.5685215 -0.1867832  0.5823712  2.0651455 

Number of Observations: 49
Number of Groups: 27 
emmeans::emmeans(Model_neutral, pairwise ~ time_point)
$emmeans
 time_point  emmean   SE df lower.CL upper.CL
 Acclimation  31.60 1.77 25    27.96     35.2
 Antibiotics   9.68 1.87 20     5.77     13.6

Results are averaged over the levels of: Population 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast                  estimate  SE df t.ratio p.value
 Acclimation - Antibiotics     21.9 2.4 20   9.147  <.0001

Results are averaged over the levels of: Population 
Degrees-of-freedom method: containment 
emmeans::emmeans(Model_neutral, pairwise ~ Population)
$emmeans
 Population emmean   SE df lower.CL upper.CL
 Cold_wet     13.4 1.61 20     10.0     16.7
 Hot_dry      27.9 2.22 20     23.3     32.5

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast           estimate   SE df t.ratio p.value
 Cold_wet - Hot_dry    -14.5 2.74 20  -5.288  <.0001

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 

Phylogenetic

Model_phylo <- nlme::lme(fixed = phylogenetic ~ time_point*Population, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_phylo)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC   logLik
  203.6546 214.4946 -95.8273

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:   0.7173285 1.688356

Fixed effects:  phylogenetic ~ time_point * Population 
                                            Value Std.Error DF   t-value p-value
(Intercept)                              5.365728 0.4446271 25 12.067930  0.0000
time_pointAntibiotics                   -0.761455 0.6038652 20 -1.260969  0.2218
PopulationHot_dry                        1.097291 0.7560383 25  1.451370  0.1591
time_pointAntibiotics:PopulationHot_dry -0.702255 1.0216421 20 -0.687379  0.4997
 Correlation: 
                                        (Intr) tm_pnA PpltH_
time_pointAntibiotics                   -0.631              
PopulationHot_dry                       -0.588  0.371       
time_pointAntibiotics:PopulationHot_dry  0.373 -0.591 -0.629

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.85993684 -0.51436620 -0.05130501  0.48944142  1.80209573 

Number of Observations: 49
Number of Groups: 27 
emmeans::emmeans(Model_neutral, pairwise ~ time_point)
$emmeans
 time_point  emmean   SE df lower.CL upper.CL
 Acclimation  31.60 1.77 25    27.96     35.2
 Antibiotics   9.68 1.87 20     5.77     13.6

Results are averaged over the levels of: Population 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast                  estimate  SE df t.ratio p.value
 Acclimation - Antibiotics     21.9 2.4 20   9.147  <.0001

Results are averaged over the levels of: Population 
Degrees-of-freedom method: containment 
emmeans::emmeans(Model_neutral, pairwise ~ Population)
$emmeans
 Population emmean   SE df lower.CL upper.CL
 Cold_wet     13.4 1.61 20     10.0     16.7
 Hot_dry      27.9 2.22 20     23.3     32.5

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast           estimate   SE df t.ratio p.value
 Cold_wet - Hot_dry    -14.5 2.74 20  -5.288  <.0001

Results are averaged over the levels of: time_point 
Degrees-of-freedom method: containment 

4.3.2 Beta diversity

samples_to_keep <- sample_metadata %>%
  filter(time_point %in% c("Acclimation", "Antibiotics")) %>% 
  select(Tube_code) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(time_point %in% c("Acclimation", "Antibiotics"))

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.03530 0.035300 10.073    999  0.004 **
Residuals 47 0.16471 0.003505                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(richness ~ time_point*Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
time_point 1 1.9348458 0.1002347 6.065418 0.001
Population 1 2.1356198 0.1106358 6.694811 0.001
time_point:Population 1 0.8778553 0.0454773 2.751930 0.004
Residual 45 14.3548318 0.7436522 NA NA
Total 48 19.3031527 1.0000000 NA NA
subset_meta_arrange <- column_to_rownames(subset_meta, "Tube_code")
subset_meta_arrange<-subset_meta_arrange[labels(richness),]

pairwise <- pairwise.adonis(richness, subset_meta_arrange$Population, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Cold_wet vs Hot_dry 1 2.135115 5.845187 0.1106096 0.001 0.001 **
pairwise <- pairwise.adonis(richness, subset_meta_arrange$time_point, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Acclimation vs Antibiotics 1 1.934846 5.235844 0.1002347 0.001 0.001 **

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.051657 0.051657 9.9224    999  0.003 **
Residuals 47 0.244689 0.005206                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(neutral ~ time_point*Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
time_point 1 2.0429592 0.11061666 7.254623 0.001
Population 1 2.8764490 0.15574623 10.214376 0.001
time_point:Population 1 0.8770558 0.04748846 3.114457 0.002
Residual 45 12.6723552 0.68614864 NA NA
Total 48 18.4688192 1.00000000 NA NA
subset_meta_arrange <- column_to_rownames(subset_meta, "Tube_code")
subset_meta_arrange<-subset_meta_arrange[labels(neutral),]

pairwise <- pairwise.adonis(neutral, subset_meta_arrange$Population, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Cold_wet vs Hot_dry 1 2.876596 8.670991 0.1557542 0.001 0.001 **
pairwise <- pairwise.adonis(neutral, subset_meta_arrange$time_point, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Acclimation vs Antibiotics 1 2.042959 5.845605 0.1106167 0.001 0.001 **

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$time_point) %>% permutest(.)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.67439 0.67439 55.932    999  0.001 ***
Residuals 47 0.56670 0.01206                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(phylo ~ time_point*Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
time_point 1 1.8783556 0.23032747 16.358708 0.001
Population 1 0.7551942 0.09260332 6.577030 0.001
time_point:Population 1 0.3545686 0.04347786 3.087959 0.022
Residual 45 5.1670341 0.63359135 NA NA
Total 48 8.1551525 1.00000000 NA NA
subset_meta_arrange <- column_to_rownames(subset_meta, "Tube_code")
subset_meta_arrange<-subset_meta_arrange[labels(phylo),]

pairwise <- pairwise.adonis(phylo, subset_meta_arrange$Population, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Cold_wet vs Hot_dry 1 0.753589 4.785297 0.09240649 0.001 0.001 **
pairwise <- pairwise.adonis(phylo, subset_meta_arrange$time_point, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Acclimation vs Antibiotics 1 1.878356 14.06493 0.2303275 0.001 0.001 **

dbRDA

#Richness
cca_ord <- capscale(formula = richness ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                 "subset_meta$time_pointAntibiotics" = "Antibiotics",
                                 "subset_meta$PopulationHot_dry"="Warm population",
                                 "subset_meta$time_pointAntibiotics:subset_meta$PopulationHot_dry"="Interaction Wam population")

beta_richness_rda <-CAP_df %>%
  group_by(Population, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
  scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()

#Neutral
cca_ord <- capscale(formula = neutral ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                  "subset_meta$time_pointAntibiotics" = "Antibiotics",
                                 "subset_meta$PopulationHot_dry"="Warm population",
                                 "subset_meta$time_pointAntibiotics:subset_meta$PopulationHot_dry"="Interaction Wam population")

beta_neutral_rda <- CAP_df %>%
  group_by(Population, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
  scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50"))+
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()


#Phylogenetic
cca_ord <- capscale(formula = phylo ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable, 
                                  "subset_meta$time_pointAntibiotics" = "Antibiotics",
                                 "subset_meta$PopulationHot_dry"="Warm population",
                                 "subset_meta$time_pointAntibiotics:subset_meta$PopulationHot_dry"="Interaction Wam population")

beta_phylogenetic_rda<- CAP_df %>%
  group_by(type, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
  scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()
ggarrange(beta_richness_rda, beta_neutral_rda, beta_phylogenetic_rda, ncol=3, nrow=1, common.legend = TRUE, legend="right")

4.3.3 Cold population

4.3.3.1 Structural zeros

acclimation_samples <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point == "Acclimation") %>% 
  dplyr::select(sample) %>%
  pull()

antibiotics_samples <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point == "Antibiotics") %>% 
  dplyr::select(sample) %>% pull()

structural_zeros <- genome_counts_filt %>% 
  select(one_of(c("genome",subset_meta$sample))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>% 
   rowwise() %>% #compute for each row (genome)
   mutate(all_zeros_acclimation = all(c_across(all_of(acclimation_samples)) == 0)) %>% 
   mutate(all_zeros_antibiotics = all(c_across(all_of(antibiotics_samples)) == 0)) %>% 
   mutate(average_acclimation = mean(c_across(all_of(acclimation_samples)), na.rm = TRUE)) %>% 
   mutate(average_antibiotics = mean(c_across(all_of(antibiotics_samples)), na.rm = TRUE)) %>% 
   filter(all_zeros_acclimation == TRUE || all_zeros_antibiotics==TRUE)  %>% 
   mutate(present = case_when(
      all_zeros_acclimation & !all_zeros_antibiotics ~ "antibiotics",
      !all_zeros_acclimation & all_zeros_antibiotics ~ "acclimation",
      !all_zeros_acclimation & !all_zeros_antibiotics ~ "None",
      TRUE ~ NA_character_
    )) %>%
   mutate(average = ifelse(present == "acclimation", average_acclimation, average_antibiotics)) %>%
   dplyr::select(genome, present, average) %>%
   left_join(genome_metadata, by=join_by(genome==genome)) %>%
   arrange(present,-average)

Acclimation

structural_zeros %>% 
  filter(present=="acclimation") %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count)) 
# A tibble: 9 × 2
# Rowwise: 
  phylum               sample_count
  <chr>                       <int>
1 p__Bacillota_A                 43
2 p__Bacteroidota                15
3 p__Bacillota                    8
4 p__Pseudomonadota               7
5 p__Cyanobacteriota              3
6 p__Verrucomicrobiota            2
7 p__Bacillota_B                  1
8 p__Bacillota_C                  1
9 p__Fusobacteriota               1

Antibiotics

structural_zeros %>% 
  filter(present=="antibiotics") %>% 
  select(1,5:10)
# A tibble: 4 × 7
# Rowwise: 
  genome                phylum            class                  order                 family                genus         species            
  <chr>                 <chr>             <chr>                  <chr>                 <chr>                 <chr>         <chr>              
1 AH1_2nd_7:bin_000003  p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales   f__Enterobacteriaceae g__Proteus    s__Proteus cibarius
2 LI1_2nd_7:bin_000001  p__Bacillota_A    c__Clostridia          o__Clostridiales      f__Clostridiaceae     g__Sarcina    s__                
3 AH1_2nd_7:bin_000055  p__Bacillota      c__Bacilli             o__Mycoplasmatales    f__Mycoplasmoidaceae  g__Ureaplasma s__                
4 AH1_2nd_13:bin_000025 p__Bacillota_A    c__Clostridia          o__Christensenellales f__UBA3700            g__           s__                

Community elements differences in structural zeros

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))
uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)
# A tibble: 0 × 3
# ℹ 3 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>

4.3.3.2 Differential abundance analysis: composition

phylo_samples <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point %in% c("Acclimation", "Antibiotics") )%>% 
  column_to_rownames("sample") %>% 
  sample_data()
phylo_genome <- genome_counts_filt %>% 
  filter(!genome %in% structural_zeros$genome) %>% 
  select(one_of(c("genome",rownames(phylo_samples)))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
  column_to_rownames("genome") %>% 
  otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>% 
  filter(genome %in% rownames(phylo_genome)) %>% 
  column_to_rownames("genome") %>% 
  dplyr::select(domain,phylum,class,order,family,genus,species) %>% 
  as.matrix() %>% 
  tax_table()

physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)
ancom_rand_output = ancombc2(data = physeq_genome_filtered, 
                  assay_name = "counts", 
                  tax_level = NULL,
                  fix_formula = "time_point",
                  rand_formula = "(1|individual)",
                  p_adj_method = "holm", 
                  pseudo_sens = TRUE,
                  prv_cut =0, 
                  lib_cut = 0, 
                  s0_perc = 0.05,
                  group = NULL, 
                  struc_zero = FALSE, 
                  neg_lb = FALSE,
                  alpha = 0.05, 
                  n_cl = 2, 
                  verbose = TRUE,
                  global = FALSE, 
                  pairwise = FALSE, 
                  dunnet = FALSE, 
                  trend = FALSE,
                  iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
                  em_control = list(tol = 1e-5, max_iter = 100),
                  lme_control = lme4::lmerControl(),
                  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100), 
                  trend_control = NULL)
save(ancom_rand_output, 
     file = "data/ancombc_antib_accl.Rdata")
load("data/ancombc_antib_accl.Rdata")
taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
  rownames_to_column(., "taxon") %>%
  mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))

ancombc_rand_table_mag <- ancom_rand_output$res %>%
  dplyr::select(taxon, lfc_time_point2_Antibiotics, p_time_point2_Antibiotics) %>%
  filter(p_time_point2_Antibiotics < 0.05) %>%
  dplyr::arrange(p_time_point2_Antibiotics) %>%
  merge(., taxonomy, by="taxon") %>%
  mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
  dplyr::arrange(lfc_time_point2_Antibiotics)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", ""))  %>%
  right_join(taxonomy, by=join_by(phylum == phylum)) %>%
  dplyr::select(phylum, colors) %>%
  mutate(colors = str_c(colors, "80"))  %>% #add 80% alpha
    unique() %>%
    dplyr::arrange(phylum)

tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
  
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
    dplyr::arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()
ancombc_rand_table_mag%>%
      mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_point2_Antibiotics, y=forcats::fct_reorder(genome,lfc_time_point2_Antibiotics), fill=phylum)) + #forcats::fct_rev()
  geom_col() + 
  scale_fill_manual(values=tax_color) + 
  geom_hline(yintercept=0) + 
#  coord_flip()+
  theme(panel.background = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 8),
        axis.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14, face = "bold"),
        legend.position = "right", legend.box = "vertical")+
  xlab("log2FoldChange") + 
  ylab("Species")+
  guides(fill=guide_legend(title="Phylum"))

Phyla of the significant MAGs

ancombc_rand_table_mag%>%
  filter(lfc_time_point2_Antibiotics<0)  %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count))  
            phylum sample_count
1     Bacteroidota           14
2      Bacillota_A            4
3        Bacillota            2
4 Campylobacterota            1
ancombc_rand_table_mag%>%
  filter(lfc_time_point2_Antibiotics<0)  %>% 
  select(phylum, genus, species)
             phylum           genus species
1  Campylobacterota  Helicobacter_J        
2       Bacillota_A                        
3      Bacteroidota     Bacteroides        
4      Bacteroidota     Odoribacter        
5      Bacteroidota     Bacteroides        
6         Bacillota   Coprobacillus        
7         Bacillota                        
8      Bacteroidota     Phocaeicola        
9      Bacteroidota     Bacteroides        
10     Bacteroidota     Phocaeicola        
11     Bacteroidota     Odoribacter        
12     Bacteroidota     Bacteroides        
13     Bacteroidota     Bacteroides        
14     Bacteroidota Parabacteroides        
15     Bacteroidota                        
16     Bacteroidota Parabacteroides        
17      Bacillota_A                        
18     Bacteroidota       Alistipes        
19      Bacillota_A                        
20      Bacillota_A                        
21     Bacteroidota     Odoribacter        
ancombc_rand_table_mag%>%
  filter(lfc_time_point2_Antibiotics<0)  %>% 
  count(genus, name = "sample_count") %>%
  arrange(desc(sample_count))
            genus sample_count
1                            6
2     Bacteroides            5
3     Odoribacter            3
4 Parabacteroides            2
5     Phocaeicola            2
6       Alistipes            1
7   Coprobacillus            1
8  Helicobacter_J            1

4.3.3.3 Differences in the functional capacity

GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>% 
  select_if(~ !is.numeric(.) || sum(.) != 0)

GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)

GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)

sample_sub <- sample_metadata %>%
  filter(Population == "Cold_wet" & treatment %in% c("Acclimation", "Antibiotics"))

genome_counts_row <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% 
    select(one_of(c("genome",sample_sub$sample))) %>%
  column_to_rownames(., "genome") 

GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)
4.3.3.3.1 Community elements differences:

MCI

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))
# A tibble: 2 × 3
  treatment     MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.335 0.0324
2 Antibiotics 0.247 0.136 
MCI <- GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.94332, p-value = 0.09304
wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 190, p-value = 0.01768
alternative hypothesis: true location shift is not equal to 0
element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))
uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

significant_elements <- element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)%>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "sample")%>% 
  left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))

difference_table <- element_gift_filt %>%
  select(-sample) %>%
  group_by(treatment) %>%
  summarise(across(everything(), mean)) %>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric) %>%
  rownames_to_column(., "Elements") %>%
  left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>% 
  arrange(Function) %>% 
  mutate(Difference=Acclimation-Antibiotics)%>% 
  mutate(group_color = ifelse(Difference <0, "Antibiotics","Acclimation")) 
difference_table %>%
  ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) + 
  geom_col() +
#  geom_point(size=4) + 
  scale_fill_manual(values=c('#008080', '#00808050')) +
  geom_hline(yintercept=0) + 
  coord_flip()+
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        legend.position = "right", 
        panel.background = element_blank(),
          panel.grid.major = element_line(size = 0.15, linetype = 'solid',
                                colour = "grey"))+
  xlab("Microbial Functional Capacity") + 
  ylab("Mean difference")+
  labs(fill="Treatment")

4.3.3.3.2 Community functions differences

MCI

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))
# A tibble: 2 × 3
  treatment     MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.330 0.0320
2 Antibiotics 0.256 0.126 
MCI <- GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.95742, p-value = 0.2331
wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 188, p-value = 0.02195
alternative hypothesis: true location shift is not equal to 0
function_gift <- GIFTs_functions_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  merge(., sample_metadata[c(12,13)], by="sample")
unique_funct_db<- GIFT_db[c(3,4,5)] %>% 
  distinct(Code_function, .keep_all = TRUE)

significant_functional <- function_gift %>%
  pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
  mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
  filter(p_adjust < 0.05)%>%
  left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))
Code_function Acclimation Antibiotics Function Difference treatment
B06 0.68110280 0.47325490 Organic anion biosynthesis 0.20784790 Acclimation
B02 0.59963040 0.41562100 Amino acid biosynthesis 0.18400940 Acclimation
D02 0.38587770 0.20636760 Polysaccharide degradation 0.17951010 Acclimation
S03 0.27103471 0.09357224 Spore 0.17746247 Acclimation
B01 0.84215140 0.69078910 Nucleic acid biosynthesis 0.15136230 Acclimation
B07 0.44558700 0.30432910 Vitamin biosynthesis 0.14125790 Acclimation
B08 0.44596150 0.32044710 Aromatic compound biosynthesis 0.12551440 Acclimation
D09 0.25533360 0.13438350 Antibiotic degradation 0.12095010 Acclimation
B04 0.54457570 0.42921780 SCFA biosynthesis 0.11535790 Acclimation
D03 0.29210750 0.20698550 Sugar degradation 0.08512200 Acclimation
D05 0.28856110 0.22258820 Amino acid degradation 0.06597290 Acclimation
B10 0.05947806 0.03621687 Antibiotic biosynthesis 0.02326119 Acclimation

4.3.4 Warm population

4.3.4.1 Structural zeros

acclimation_samples <- sample_metadata %>%
  filter(Population == "Hot_dry" & time_point == "Acclimation") %>% 
  dplyr::select(sample) %>%
  pull()

antibiotics_samples <- sample_metadata %>%
  filter(Population == "Hot_dry" & time_point == "Antibiotics") %>% 
  dplyr::select(sample) %>% pull()

structural_zeros <- genome_counts_filt %>% 
  select(one_of(c("genome",subset_meta$sample))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>% 
   rowwise() %>% #compute for each row (genome)
   mutate(all_zeros_acclimation = all(c_across(all_of(acclimation_samples)) == 0)) %>% 
   mutate(all_zeros_antibiotics = all(c_across(all_of(antibiotics_samples)) == 0)) %>% 
   mutate(average_acclimation = mean(c_across(all_of(acclimation_samples)), na.rm = TRUE)) %>% 
   mutate(average_antibiotics = mean(c_across(all_of(antibiotics_samples)), na.rm = TRUE)) %>% 
   filter(all_zeros_acclimation == TRUE || all_zeros_antibiotics==TRUE)  %>% 
   mutate(present = case_when(
      all_zeros_acclimation & !all_zeros_antibiotics ~ "antibiotics",
      !all_zeros_acclimation & all_zeros_antibiotics ~ "acclimation",
      !all_zeros_acclimation & !all_zeros_antibiotics ~ "None",
      TRUE ~ NA_character_
    )) %>%
   mutate(average = ifelse(present == "acclimation", average_acclimation, average_antibiotics)) %>%
   dplyr::select(genome, present, average) %>%
   left_join(genome_metadata, by=join_by(genome==genome)) %>%
   arrange(present,-average)

Acclimation

structural_zeros %>% 
  filter(present=="acclimation") %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count)) 
# A tibble: 12 × 2
# Rowwise: 
   phylum               sample_count
   <chr>                       <int>
 1 p__Bacillota_A                 63
 2 p__Bacteroidota                32
 3 p__Bacillota                   16
 4 p__Cyanobacteriota              9
 5 p__Pseudomonadota               7
 6 p__Bacillota_B                  3
 7 p__Desulfobacterota             3
 8 p__Bacillota_C                  2
 9 p__Verrucomicrobiota            2
10 p__Actinomycetota               1
11 p__Campylobacterota             1
12 p__Elusimicrobiota              1

Antibiotics

structural_zeros %>% 
  filter(present=="antibiotics") %>% 
  select(1,5:10)
# A tibble: 7 × 7
# Rowwise: 
  genome                phylum            class                  order                 family                genus              species                       
  <chr>                 <chr>             <chr>                  <chr>                 <chr>                 <chr>              <chr>                         
1 AH1_2nd_7:bin_000003  p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales   f__Enterobacteriaceae g__Proteus         s__Proteus cibarius           
2 LI1_2nd_8:bin_000030  p__Actinomycetota c__Actinomycetia       o__Mycobacteriales    f__Mycobacteriaceae   g__Corynebacterium s__                           
3 LI1_2nd_8:bin_000077  p__Bacteroidota   c__Bacteroidia         o__Bacteroidales      f__Tannerellaceae     g__Parabacteroides s__Parabacteroides sp003480915
4 LI1_2nd_5:bin_000032  p__Bacillota_A    c__Clostridia          o__Christensenellales f__UBA1242            g__Caccovivens     s__                           
5 AH1_2nd_18:bin_000013 p__Bacteroidota   c__Bacteroidia         o__Bacteroidales      f__Tannerellaceae     g__Parabacteroides s__                           
6 LI1_2nd_5:bin_000069  p__Bacillota      c__Bacilli             o__RF39               f__UBA660             g__                s__                           
7 AH1_2nd_20:bin_000061 p__Bacillota      c__Bacilli             o__Lactobacillales    f__Enterococcaceae    g__Enterococcus    s__                           

Community elements differences in structural zeros

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))
uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)
# A tibble: 65 × 3
   trait p_value p_adjust
   <chr>   <dbl>    <dbl>
 1 B0101 0.00280   0.0184
 2 B0103 0.0196    0.0437
 3 B0105 0.00559   0.0184
 4 B0106 0.00559   0.0184
 5 B0205 0.00280   0.0184
 6 B0207 0.0112    0.0280
 7 B0212 0.00559   0.0184
 8 B0214 0.00662   0.0184
 9 B0216 0.00559   0.0184
10 B0307 0.0196    0.0437
# ℹ 55 more rows

4.3.4.2 Differential abundance analysis: composition

phylo_samples <- sample_metadata %>%
  filter(Population == "Hot_dry" & time_point %in% c("Acclimation", "Antibiotics") )%>% 
  column_to_rownames("sample") %>% 
  sample_data()
phylo_genome <- genome_counts_filt %>% 
  filter(!genome %in% structural_zeros$genome) %>% 
  select(one_of(c("genome",rownames(phylo_samples)))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
  column_to_rownames("genome") %>% 
  otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>% 
  filter(genome %in% rownames(phylo_genome)) %>% 
  column_to_rownames("genome") %>% 
  dplyr::select(domain,phylum,class,order,family,genus,species) %>% 
  as.matrix() %>% 
  tax_table()

physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)
ancom_rand_output = ancombc2(data = physeq_genome_filtered, 
                  assay_name = "counts", 
                  tax_level = NULL,
                  fix_formula = "time_point",
                  rand_formula = "(1|individual)",
                  p_adj_method = "holm", 
                  pseudo_sens = TRUE,
                  prv_cut =0, 
                  lib_cut = 0, 
                  s0_perc = 0.05,
                  group = NULL, 
                  struc_zero = FALSE, 
                  neg_lb = FALSE,
                  alpha = 0.05, 
                  n_cl = 2, 
                  verbose = TRUE,
                  global = FALSE, 
                  pairwise = FALSE, 
                  dunnet = FALSE, 
                  trend = FALSE,
                  iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
                  em_control = list(tol = 1e-5, max_iter = 100),
                  lme_control = lme4::lmerControl(),
                  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100), 
                  trend_control = NULL)
save(ancom_rand_output, 
     file = "data/ancombc_antib_accl_warm.Rdata")
load("data/ancombc_antib_accl_warm.Rdata")
taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
  rownames_to_column(., "taxon") %>%
  mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))

ancombc_rand_table_mag <- ancom_rand_output$res %>%
  dplyr::select(taxon, lfc_time_point2_Antibiotics, p_time_point2_Antibiotics) %>%
  filter(p_time_point2_Antibiotics < 0.05) %>%
  dplyr::arrange(p_time_point2_Antibiotics) %>%
  merge(., taxonomy, by="taxon") %>%
  mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
  dplyr::arrange(lfc_time_point2_Antibiotics)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", ""))  %>%
  right_join(taxonomy, by=join_by(phylum == phylum)) %>%
  dplyr::select(phylum, colors) %>%
  mutate(colors = str_c(colors, "80"))  %>% #add 80% alpha
    unique() %>%
    dplyr::arrange(phylum)

tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
  
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
    dplyr::arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()
ancombc_rand_table_mag%>%
      mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_point2_Antibiotics, y=forcats::fct_reorder(genome,lfc_time_point2_Antibiotics), fill=phylum)) + #forcats::fct_rev()
  geom_col() + 
  scale_fill_manual(values=tax_color) + 
  geom_hline(yintercept=0) + 
#  coord_flip()+
  theme(panel.background = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 8),
        axis.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14, face = "bold"),
        legend.position = "right", legend.box = "vertical")+
  xlab("log2FoldChange") + 
  ylab("Species")+
  guides(fill=guide_legend(title="Phylum"))

Phyla of the significant MAGs

ancombc_rand_table_mag%>%
  filter(lfc_time_point2_Antibiotics<0)  %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count))  

ancombc_rand_table_mag%>%
  filter(lfc_time_point2_Antibiotics<0)  %>% 
  select(phylum, genus, species)

ancombc_rand_table_mag%>%
  filter(lfc_time_point2_Antibiotics<0)  %>% 
  count(genus, name = "sample_count") %>%
  arrange(desc(sample_count))

4.3.4.3 Differences in the functional capacity

GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>% 
  select_if(~ !is.numeric(.) || sum(.) != 0)

GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)

GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)

sample_sub <- sample_metadata %>%
  filter(Population == "Cold_wet" & treatment %in% c("Acclimation", "Antibiotics"))

genome_counts_row <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% 
    select(one_of(c("genome",sample_sub$sample))) %>%
  column_to_rownames(., "genome") 

GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)
4.3.4.3.1 Community elements differences:

MCI

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))
# A tibble: 2 × 3
  treatment     MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.335 0.0324
2 Antibiotics 0.247 0.136 
MCI <- GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.94332, p-value = 0.09304
wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 190, p-value = 0.01768
alternative hypothesis: true location shift is not equal to 0
element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))
uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

significant_elements <- element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)%>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "sample")%>% 
  left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))

difference_table <- element_gift_filt %>%
  select(-sample) %>%
  group_by(treatment) %>%
  summarise(across(everything(), mean)) %>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric) %>%
  rownames_to_column(., "Elements") %>%
  left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>% 
  arrange(Function) %>% 
  mutate(Difference=Acclimation-Antibiotics)%>% 
  mutate(group_color = ifelse(Difference <0, "Antibiotics","Acclimation")) 
difference_table %>%
  ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) + 
  geom_col() +
#  geom_point(size=4) + 
  scale_fill_manual(values=c("#d57d2c50","#d57d2c")) +
  geom_hline(yintercept=0) + 
  coord_flip()+
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        legend.position = "right", 
        panel.background = element_blank(),
          panel.grid.major = element_line(size = 0.15, linetype = 'solid',
                                colour = "grey"))+
  xlab("Microbial Functional Capacity") + 
  ylab("Mean difference")+
  labs(fill="Treatment")

4.3.4.3.2 Community functions differences

MCI

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))
# A tibble: 2 × 3
  treatment     MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.330 0.0320
2 Antibiotics 0.256 0.126 
MCI <- GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.95742, p-value = 0.2331
wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 188, p-value = 0.02195
alternative hypothesis: true location shift is not equal to 0
function_gift <- GIFTs_functions_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  merge(., sample_metadata[c(12,13)], by="sample")
unique_funct_db<- GIFT_db[c(3,4,5)] %>% 
  distinct(Code_function, .keep_all = TRUE)

significant_functional <- function_gift %>%
  pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
  mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
  filter(p_adjust < 0.05)%>%
  left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))
Code_function Acclimation Antibiotics Function Difference treatment
B06 0.68110280 0.47325490 Organic anion biosynthesis 0.20784790 Acclimation
B02 0.59963040 0.41562100 Amino acid biosynthesis 0.18400940 Acclimation
D02 0.38587770 0.20636760 Polysaccharide degradation 0.17951010 Acclimation
S03 0.27103471 0.09357224 Spore 0.17746247 Acclimation
B01 0.84215140 0.69078910 Nucleic acid biosynthesis 0.15136230 Acclimation
B07 0.44558700 0.30432910 Vitamin biosynthesis 0.14125790 Acclimation
B08 0.44596150 0.32044710 Aromatic compound biosynthesis 0.12551440 Acclimation
D09 0.25533360 0.13438350 Antibiotic degradation 0.12095010 Acclimation
B04 0.54457570 0.42921780 SCFA biosynthesis 0.11535790 Acclimation
D03 0.29210750 0.20698550 Sugar degradation 0.08512200 Acclimation
D05 0.28856110 0.22258820 Amino acid degradation 0.06597290 Acclimation
B10 0.05947806 0.03621687 Antibiotic biosynthesis 0.02326119 Acclimation

4.4 Does antibiotic administration remove the differences found in the two populations?

4.4.1 Shapiro test

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point == "Antibiotics" ) %>%
  filter(metric=="richness") %>% 
  summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
  pull(shapiro_p_value)
[1] 0.001165318
alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point == "Antibiotics" ) %>%
  filter(metric=="neutral") %>% 
  summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
  pull(shapiro_p_value)
[1] 0.0003462674
alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point == "Antibiotics" ) %>%
  filter(metric=="phylogenetic") %>% 
  summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
  pull(shapiro_p_value)
[1] 0.06628442

4.4.2 Alpha diversity

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  filter(time_point == "Antibiotics" )
alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness", "neutral"))) %>%
  filter(metric!="phylogenetic") %>%
  filter(time_point == "Antibiotics" ) %>% 
  ggplot(aes(y = value, x = Population, color=Population, fill=Population)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(alpha=0.5) +
  scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
  stat_compare_means(method = "t.test", show.legend = F, size = 3, label.x = c(1.5))+
  facet_grid( ~ metric, scales = "free_y")+
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral"))) %>%
  filter(metric!="phylogenetic") %>%
  filter(time_point == "Antibiotics" ) %>% 
  ggplot(aes(y = value, x = Population, color=Population, fill=Population)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(alpha=0.5) +
  scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
  stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

4.4.3 Beta diversity

samples_to_keep <- sample_metadata %>%
  filter(time_point == "Antibiotics") %>% 
  select(Tube_code) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(time_point == "Antibiotics")

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$Population) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.015377 0.0153774 6.8934    999  0.021 *
Residuals 21 0.046845 0.0022307                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(richness ~ Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 1.361743 0.1533845 3.80465 0.001
Residual 21 7.516224 0.8466155 NA NA
Total 22 8.877966 1.0000000 NA NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$Population) %>% permutest(.)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.030649 0.0306488 3.8694    999   0.06 .
Residuals 21 0.166339 0.0079209                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(neutral ~ Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 1.787698 0.208772 5.541022 0.001
Residual 21 6.775221 0.791228 NA NA
Total 22 8.562919 1.000000 NA NA

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$Population) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.012165 0.012165 0.9979    999   0.33
Residuals 21 0.256012 0.012191                     
adonis2(phylo ~ Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 0.8970751 0.1890846 4.896659 0.001
Residual 21 3.8472307 0.8109154 NA NA
Total 22 4.7443057 1.0000000 NA NA

4.5 Are the microbial communities similar in both donor samples?

4.5.1 Alpha diversity

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  filter(type =="Hot_control" & time_point %in% c( "Transplant1", "Transplant2"))
samples_to_keep <- sample_metadata %>%
  filter(type =="Hot_control" & time_point %in% c( "Transplant1", "Transplant2"))%>% 
  select(sample) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(type =="Hot_control" & time_point %in% c( "Transplant1", "Transplant2"))
alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
  filter(type =="Hot_control" & time_point %in% c( "Transplant1", "Transplant2")) %>% 
  ggplot(aes(y = value, x = time_point, color=time_point, fill=time_point)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(alpha=0.5) +
  scale_color_manual(name="time_point",
          breaks=c("Transplant1","Transplant2"),
          values=c('#d5992c', "#d5b52c")) +
      scale_fill_manual(name="time_point",
          breaks=c("Transplant1","Transplant2"),
          values=c('#d5992c50', "#d5b52c50")) +
  facet_wrap(. ~ metric, scales = "free") +
  theme_classic() +
  theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

Richness

Modelq0GLMMNB <- glmer.nb(richness ~ time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(17.9668)  ( log )
Formula: richness ~ time_point + (1 | individual)
   Data: alpha_div_meta

     AIC      BIC   logLik deviance df.resid 
   150.4    153.3    -71.2    142.4       11 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.02956 -0.52933  0.03513  0.56378  1.56130 

Random effects:
 Groups     Name        Variance Std.Dev.
 individual (Intercept) 0.009989 0.09995 
Number of obs: 15, groups:  individual, 8

Fixed effects:
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)            4.60703    0.09905  46.514   <2e-16 ***
time_pointTransplant2  0.05482    0.13299   0.412     0.68    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr)
tm_pntTrns2 -0.624

Neutral

Model_neutral <- nlme::lme(fixed = neutral ~ time_point, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_neutral)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  115.9183 118.1781 -53.95914

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    9.533528 10.15744

Fixed effects:  neutral ~ time_point 
                         Value Std.Error DF  t-value p-value
(Intercept)           43.94053  4.925212  7 8.921551  0.0000
time_pointTransplant2  4.31623  5.338413  6 0.808523  0.4497
 Correlation: 
                      (Intr)
time_pointTransplant2 -0.491

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.57714267 -0.56764018  0.05406399  0.55332481  1.12750552 

Number of Observations: 15
Number of Groups: 8 

Phylogenetic

Model_phylo <-  nlme::lme(fixed = phylogenetic ~ time_point, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_phylo)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  41.16537 43.42517 -16.58269

Random effects:
 Formula: ~1 | individual
        (Intercept)  Residual
StdDev:     1.16159 0.3042944

Fixed effects:  phylogenetic ~ time_point 
                          Value Std.Error DF   t-value p-value
(Intercept)            5.785003 0.4245418  7 13.626462  0.0000
time_pointTransplant2 -0.018098 0.1623255  6 -0.111494  0.9149
 Correlation: 
                      (Intr)
time_pointTransplant2 -0.168

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.2428123 -0.3756145 -0.0502414  0.4569146  1.3175579 

Number of Observations: 15
Number of Groups: 8 

4.5.2 Beta diversity

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$treatment) %>% permutest(.)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00014 0.0001396 0.0173    999  0.907
Residuals 13 0.10481 0.0080621                     
adonis2(richness ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(richness))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 0.06889268 0.03010712 0.4035421 0.9375
Residual 13 2.21935895 0.96989288 NA NA
Total 14 2.28825162 1.00000000 NA NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000655 0.0006546 0.0514    999  0.835
Residuals 13 0.165409 0.0127237                     
adonis2(neutral ~ treatment,
        data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 0.09294671 0.03882177 0.5250671 0.7265625
Residual 13 2.30124351 0.96117823 NA NA
Total 14 2.39419022 1.00000000 NA NA

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.003691 0.0036912 0.3071    999  0.674
Residuals 13 0.156266 0.0120205                     
adonis2(phylo ~ treatment,
        data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 0.009773632 0.01921019 0.2546239 0.7734375
Residual 13 0.498999565 0.98078981 NA NA
Total 14 0.508773196 1.00000000 NA NA

4.6 Does the donor sample maintain the microbial community found in acclimation?

sample_metadata <- sample_metadata %>%
    mutate(treatment=case_when(
    treatment %in% c("Transplant1", "Transplant2") ~ "Transplant",
    TRUE ~ treatment
  ))

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant"))
sample_metadata <- sample_metadata %>%
    mutate(treatment=case_when(
    treatment %in% c("Transplant1", "Transplant2") ~ "Transplant",
    TRUE ~ treatment
  ))
samples_to_keep <- sample_metadata %>%
  filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant")) %>% 
  select(sample) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant"))

4.6.1 Alpha diversity

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
  filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant")) %>% 
  ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(alpha=0.5) +
  scale_color_manual(name="treatment",
          breaks=c("Acclimation","Transplant"),
          values=c("#d57d2c", "#d5b52c")) +
      scale_fill_manual(name="treatment",
          breaks=c("Acclimation","Transplant"),
          values=c("#d57d2c50", "#d5b52c50")) +
  facet_wrap(. ~ metric, scales = "free") +
  theme_classic() +
  theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

Richness

Modelq0GLMMNB <- glmer.nb(richness ~ treatment+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(19.6072)  ( log )
Formula: richness ~ treatment + (1 | individual)
   Data: alpha_div_meta

     AIC      BIC   logLik deviance df.resid 
   228.0    232.7   -110.0    220.0       20 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.37872 -0.53180  0.06467  0.46904  1.76837 

Random effects:
 Groups     Name        Variance Std.Dev.
 individual (Intercept) 0        0       
Number of obs: 24, groups:  individual, 9

Fixed effects:
                    Estimate Std. Error z value Pr(>|z|)    
(Intercept)           4.4569     0.0834  53.441   <2e-16 ***
treatmentTransplant   0.1855     0.1049   1.769   0.0769 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr)
trtmntTrnsp -0.795
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

Neutral

Model_neutral <- nlme::lme(fixed = neutral ~ treatment, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_neutral)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  185.3367 189.7009 -88.66837

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev: 0.001909967 12.18199

Fixed effects:  neutral ~ treatment 
                       Value Std.Error DF   t-value p-value
(Intercept)         44.09058  4.060663 14 10.857976  0.0000
treatmentTransplant  1.80350  5.136377 14  0.351122  0.7307
 Correlation: 
                    (Intr)
treatmentTransplant -0.791

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.66610957 -0.48314823 -0.03902931  0.62479437  2.02597977 

Number of Observations: 24
Number of Groups: 9 

Phylogenetic

Model_phylo <- nlme::lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_phylo)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
      AIC      BIC   logLik
  80.6646 85.02877 -36.3323

Random effects:
 Formula: ~1 | individual
        (Intercept)  Residual
StdDev:   0.7197786 0.9593315

Fixed effects:  phylogenetic ~ treatment 
                        Value Std.Error DF   t-value p-value
(Intercept)          6.463020 0.3997775 14 16.166544  0.0000
treatmentTransplant -0.577492 0.4112749 14 -1.404151  0.1821
 Correlation: 
                    (Intr)
treatmentTransplant -0.622

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max 
-2.329685717 -0.568101055 -0.001989187  0.552370170  1.529548991 

Number of Observations: 24
Number of Groups: 9 

4.6.2 Beta diversity

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00836 0.0083600 1.4409    999  0.273
Residuals 22 0.12764 0.0058019                     
adonis2(richness ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(richness))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 0.2657351 0.06396199 1.503319 0.111
Residual 22 3.8888430 0.93603801 NA NA
Total 23 4.1545781 1.00000000 NA NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000661 0.0006614 0.0736    999  0.798
Residuals 22 0.197804 0.0089911                     
adonis2(neutral ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(neutral))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 0.3164816 0.07596593 1.808646 0.078
Residual 22 3.8496178 0.92403407 NA NA
Total 23 4.1660995 1.00000000 NA NA

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00204 0.0020405 0.2622    999  0.683
Residuals 22 0.17118 0.0077807                     
adonis2(phylo ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(phylo))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 0.0412056 0.056244 1.31111 0.235
Residual 22 0.6914166 0.943756 NA NA
Total 23 0.7326222 1.000000 NA NA
p0<-richness %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 4) +
  scale_color_manual(name="treatment",
          breaks=c("Acclimation","Transplant"),
          values=c("#d57d2c", "#d5b52c")) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )

p1<-neutral %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 4) +
  scale_color_manual(name="treatment",
          breaks=c("Acclimation","Transplant"),
          values=c("#d57d2c", "#d5b52c")) +  
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )

p2<-phylo %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 4) +
  scale_color_manual(name="treatment",
          breaks=c("Acclimation","Transplant"),
          values=c("#d57d2c", "#d5b52c")) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )
ggarrange(p0, p1, p2, ncol=3, nrow=1, common.legend = TRUE, legend="right")

4.7 Is the donor sample microbiota different to recipient microbiota?

samples_to_keep <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("Acclimation", "Transplant")) %>% 
  select(sample) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("Acclimation", "Transplant"))

4.7.1 Alpha diversity

sample_metadata <- sample_metadata %>%
    mutate(treatment=case_when(
    treatment %in% c("Transplant1", "Transplant2") ~ "Transplant",
    TRUE ~ treatment
  ))

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  filter(type == "Treatment" & treatment %in% c("Acclimation","Transplant"))
alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
  filter(type == "Treatment" & treatment %in% c("Acclimation", "Transplant")) %>% 
  ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(alpha=0.5) +
  scale_color_manual(name="treatment",
          breaks=c("Acclimation","Transplant"),
          values=c('#008080', "#d5b52c")) +
      scale_fill_manual(name="treatment",
          breaks=c("Acclimation","Transplant"),
          values=c('#00808050', "#d5b52c50")) +
  facet_wrap(. ~ metric, scales = "free") +
  theme_classic() +
  theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

Richness

Modelq0GLMMNB <- glmer.nb(richness ~ treatment+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(6.4835)  ( log )
Formula: richness ~ treatment + (1 | individual)
   Data: alpha_div_meta

     AIC      BIC   logLik deviance df.resid 
   219.3    223.6   -105.6    211.3       18 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.9573 -0.5447  0.1304  0.5773  1.2167 

Random effects:
 Groups     Name        Variance  Std.Dev. 
 individual (Intercept) 5.302e-13 7.281e-07
Number of obs: 22, groups:  individual, 9

Fixed effects:
                    Estimate Std. Error z value Pr(>|z|)    
(Intercept)           3.8067     0.1400   27.19  < 2e-16 ***
treatmentTransplant   0.8399     0.1795    4.68 2.87e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr)
trtmntTrnsp -0.780
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

Neutral

Model_neutral <-nlme:: lme(fixed = neutral ~ treatment, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_neutral)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
      AIC      BIC    logLik
  169.583 173.5659 -80.79148

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    4.713488 11.40939

Fixed effects:  neutral ~ treatment 
                       Value Std.Error DF  t-value p-value
(Intercept)         17.31015  4.114893 12 4.206708  0.0012
treatmentTransplant 28.70744  5.016259 12 5.722879  0.0001
 Correlation: 
                    (Intr)
treatmentTransplant -0.701

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-2.4861810 -0.5958348  0.1169620  0.5490436  1.7889456 

Number of Observations: 22
Number of Groups: 9 

Phylogenetic

Model_phylo <- nlme::lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_phylo)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  80.93175 84.91468 -36.46587

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:   0.8689563 1.118803

Fixed effects:  phylogenetic ~ treatment 
                       Value Std.Error DF   t-value p-value
(Intercept)         5.488256 0.4722057 12 11.622595    0.00
treatmentTransplant 0.226402 0.5020134 12  0.450988    0.66
 Correlation: 
                    (Intr)
treatmentTransplant -0.587

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.42844864 -0.44112019 -0.05260838  0.51238319  1.55109109 

Number of Observations: 22
Number of Groups: 9 

4.7.2 Beta diversity

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.12795 0.127946 13.179    999  0.001 ***
Residuals 20 0.19416 0.009708                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(richness ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(richness))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 1.777815 0.2760196 7.625059 0.001
Residual 20 4.663086 0.7239804 NA NA
Total 21 6.440902 1.0000000 NA NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.071502 0.071502 5.4967    999  0.029 *
Residuals 20 0.260166 0.013008                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(neutral ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(neutral))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 2.112572 0.3210813 9.458609 0.001
Residual 20 4.466983 0.6789187 NA NA
Total 21 6.579556 1.0000000 NA NA

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.04731 0.047307 2.6157    999  0.117
Residuals 20 0.36172 0.018086                     
adonis2(phylo ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(phylo))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 0.3778248 0.239656 6.303884 0.001
Residual 20 1.1987049 0.760344 NA NA
Total 21 1.5765298 1.000000 NA NA
p0<-richness %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 4) +
  scale_color_manual(name="treatment",
          breaks=c("Acclimation","Transplant"),
          values=c('#008080', "#d5b52c")) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )

p1<-neutral %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 4) +
  scale_color_manual(name="treatment",
          breaks=c("Acclimation","Transplant"),
          values=c('#008080', "#d5b52c")) +  
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )

p2<-phylo %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 4) +
  scale_color_manual(name="treatment",
          breaks=c("Acclimation","Transplant"),
          values=c('#008080', "#d5b52c")) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )
ggarrange(p0, p1, p2, ncol=3, nrow=1, common.legend = TRUE, legend="right")

4.8 Does FMT change the microbial community over time?

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  filter(type == "Treatment" & treatment %in% c("FMT1","FMT2") ) 

4.8.1 Alpha diversity

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral", "phylogenetic"))) %>%
  filter(type=="Treatment" & treatment %in% c("FMT1", "FMT2")) %>% 
  ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(alpha=0.5) +
  scale_color_manual(name="treatment",
          values=c("#76b183", "#2b4b32")) +
      scale_fill_manual(name="treatment",
          values=c("#76b18350", "#2b4b3250")) +
  facet_wrap(. ~ metric, scales = "free") +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

Richness

Modelq0GLMMNB <- glmer.nb(richness ~ treatment+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(4.3876)  ( log )
Formula: richness ~ treatment + (1 | individual)
   Data: alpha_div_meta

     AIC      BIC   logLik deviance df.resid 
   171.0    174.3    -81.5    163.0       13 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.73495 -0.71265 -0.07086  0.40744  1.88240 

Random effects:
 Groups     Name        Variance  Std.Dev. 
 individual (Intercept) 1.073e-11 3.275e-06
Number of obs: 17, groups:  individual, 9

Fixed effects:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)     3.9343     0.1759  22.369   <2e-16 ***
treatmentFMT2   0.4052     0.2402   1.687   0.0917 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr)
tretmntFMT2 -0.732
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

Neutral

Model_neutral <- nlme::lme(fixed = neutral ~ treatment, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_neutral)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC   logLik
  128.7946 131.6268 -60.3973

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    3.497472 11.25384

Fixed effects:  neutral ~ treatment 
                 Value Std.Error DF  t-value p-value
(Intercept)   24.65947  4.164756  8 5.920988  0.0004
treatmentFMT2 14.87494  5.482531  7 2.713151  0.0301
 Correlation: 
              (Intr)
treatmentFMT2 -0.7  

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.76462324 -0.55701822 -0.04763166  0.55267588  1.30333436 

Number of Observations: 17
Number of Groups: 9 

Phylogenetic

Model_phylo <- nlme::lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_phylo)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC   logLik
  56.05281 58.88501 -24.0264

Random effects:
 Formula: ~1 | individual
        (Intercept)  Residual
StdDev:   0.6954001 0.8350046

Fixed effects:  phylogenetic ~ treatment 
                 Value Std.Error DF   t-value p-value
(Intercept)   4.154090 0.3805932  8 10.914778  0.0000
treatmentFMT2 0.912624 0.4105974  7  2.222673  0.0616
 Correlation: 
              (Intr)
treatmentFMT2 -0.583

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.1361565 -0.5779234 -0.1465750  0.4815529  2.3352020 

Number of Observations: 17
Number of Groups: 9 

4.8.2 Beta diversity

samples_to_keep <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("FMT1", "FMT2")) %>% 
  select(sample) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("FMT1", "FMT2"))

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.017610 0.017610 2.8225    999  0.115
Residuals 15 0.093585 0.006239                     
adonis2(richness ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(richness))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 0.3560084 0.07968734 1.298809 0.00390625
Residual 15 4.1115570 0.92031266 NA NA
Total 16 4.4675655 1.00000000 NA NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.009828 0.0098278 0.7113    999  0.442
Residuals 15 0.207263 0.0138175                     
adonis2(neutral ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(neutral))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 0.3149954 0.08110541 1.323962 0.00390625
Residual 15 3.5687823 0.91889459 NA NA
Total 16 3.8837776 1.00000000 NA NA

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.010955 0.010955 0.6602    999   0.52
Residuals 15 0.248892 0.016593                     
adonis2(phylo ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(phylo))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 0.06391536 0.09449338 1.565312 0.0703125
Residual 15 0.61248501 0.90550662 NA NA
Total 16 0.67640037 1.00000000 NA NA
p0<-richness %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 4) +
  scale_color_manual(name="treatment",
          breaks=c("FMT1", "FMT2"),
          values=c("#76b183", "#2b4b32")) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )

p1<-neutral %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 4) +
  scale_color_manual(name="treatment",
          breaks=c("FMT1", "FMT2"),
          values=c("#76b183", "#2b4b32")) +  
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )

p2<-phylo %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 4) +
  scale_color_manual(name="treatment",
          breaks=c("FMT1", "FMT2"),
          values=c("#76b183", "#2b4b32")) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )
ggarrange(p0, p1, p2, ncol=3, nrow=1, common.legend = TRUE, legend="right")

4.9 Do FMT change the microbial community compared to antibiotics and acclimation?

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  filter(type == "Treatment" & treatment %in% c("Antibiotics", "Acclimation","FMT1") ) 

4.9.1 Alpha diversity

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral", "phylogenetic"))) %>%
  filter(type=="Treatment" & treatment %in% c( "Antibiotics","Acclimation", "FMT1")) %>% 
  ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
  geom_jitter(width = 0.2, show.legend = FALSE) +
  geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
  scale_color_manual(name="Treatment",
          breaks=c("Antibiotics","Acclimation", "FMT1"),
          values=c("#003a45",'#008080', "#76b183")) +
      scale_fill_manual(name="Treatment",
          breaks=c("Antibiotics","Acclimation", "FMT1"),
          values=c("#003a4550",'#00808050',"#76b18350")) +
  facet_wrap(. ~ metric, scales = "free") +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

Richness: Problems to calculate

Modelq0GLMMNB <- glmer.nb(richness ~ treatment+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)

Neutral

Model_neutral <- lme(fixed = neutral ~ treatment, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_neutral)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC     BIC   logLik
  173.0404 178.263 -81.5202

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    4.263247 9.322271

Fixed effects:  neutral ~ treatment 
                         Value Std.Error DF   t-value p-value
(Intercept)          17.310151  3.416951 13  5.065963  0.0002
treatmentAntibiotics -8.868175  4.744329 13 -1.869216  0.0843
treatmentFMT1         7.289358  4.552796 13  1.601073  0.1334
 Correlation: 
                     (Intr) trtmnA
treatmentAntibiotics -0.596       
treatmentFMT1        -0.621  0.457

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.93185401 -0.56384573 -0.04015247  0.47657219  1.55593253 

Number of Observations: 24
Number of Groups: 9 

Phylogenetic

Model_phylo <- lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_phylo)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
      AIC      BIC   logLik
  92.7946 98.01721 -41.3973

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:   0.5507824 1.405346

Fixed effects:  phylogenetic ~ treatment 
                         Value Std.Error DF   t-value p-value
(Intercept)           5.488256 0.5031412 13 10.907984  0.0000
treatmentAntibiotics -1.198294 0.7137109 13 -1.678962  0.1170
treatmentFMT1        -1.351708 0.6855445 13 -1.971729  0.0703
 Correlation: 
                     (Intr) trtmnA
treatmentAntibiotics -0.611       
treatmentFMT1        -0.636  0.456

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.00120146 -0.44426020  0.00785757  0.38719755  1.97397308 

Number of Observations: 24
Number of Groups: 9 

4.9.2 Beta diversity

samples_to_keep <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("Acclimation", "Antibiotics", "FMT1")) %>% 
  select(sample) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("Acclimation", "Antibiotics", "FMT1"))

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq    F N.Perm Pr(>F)
Groups     2 0.01582 0.0079101 1.05    999  0.397
Residuals 21 0.15821 0.0075336                   
adonis2(richness ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(richness))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 2 1.840841 0.2018508 2.655435 0.001
Residual 21 7.278967 0.7981492 NA NA
Total 23 9.119808 1.0000000 NA NA
pairwise.adonis(richness, subset_meta$treatment, perm = 999)
                       pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics  1 0.8287871 2.293722 0.1407734   0.003      0.009   *
2        Acclimation vs FMT1  1 0.9572484 2.939042 0.1638350   0.001      0.003   *
3        Antibiotics vs FMT1  1 0.9764237 2.751191 0.1746656   0.006      0.018   .

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$treatment) %>% permutest(.)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.01796 0.0089778 0.5959    999  0.537
Residuals 21 0.31640 0.0150667                     
adonis2(neutral ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(neutral))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 2 2.357142 0.2699241 3.882065 0.001
Residual 21 6.375470 0.7300759 NA NA
Total 23 8.732613 1.0000000 NA NA
pairwise.adonis(neutral, subset_meta$treatment, perm = 999)
                       pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics  1 0.8814932 2.747044 0.1640316   0.001      0.003   *
2        Acclimation vs FMT1  1 1.3133104 4.634354 0.2360329   0.001      0.003   *
3        Antibiotics vs FMT1  1 1.3427497 4.355528 0.2509591   0.001      0.003   *

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     2 0.15410 0.077052 2.6434    999  0.102
Residuals 21 0.61214 0.029149                     
adonis2(phylo ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(phylo))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 2 1.327076 0.3647564 6.029092 0.001
Residual 21 2.311178 0.6352436 NA NA
Total 23 3.638254 1.0000000 NA NA
pairwise.adonis(phylo, subset_meta$treatment, perm = 999)
                       pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics  1 0.6885926 5.124832 0.2679674   0.004      0.012   .
2        Acclimation vs FMT1  1 0.2548027 3.292748 0.1800029   0.040      0.120    
3        Antibiotics vs FMT1  1 1.1000473 9.048069 0.4103792   0.004      0.012   .
p0<-richness %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 4) +
  scale_color_manual(name="treatment",
          breaks=c("Acclimation", "Antibiotics", "FMT1"),
          values=c("#003a45",'#008080', "#76b183")) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )

p1<-neutral %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 4) +
  scale_color_manual(name="treatment",
          breaks=c("Acclimation", "Antibiotics", "FMT1"),
          values=c("#003a45",'#008080', "#76b183")) +  
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )

p2<-phylo %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 4) +
  scale_color_manual(name="treatment",
          breaks=c("Acclimation", "Antibiotics", "FMT1"),
          values=c("#003a45",'#008080', "#76b183")) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )
ggarrange(p0, p1, p2, ncol=3, nrow=1, common.legend = TRUE, legend="right")

4.10 Is the gut microbiota similar to the donor after FMT?

sample_metadata <- sample_metadata %>%
    mutate(treatment=case_when(
    treatment %in% c("FMT1", "FMT2") ~ "FMT",
    TRUE ~ treatment
  ))

samples_to_keep <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("FMT", "Transplant")) %>% 
  select(sample) %>% 
  pull()

subset_meta <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("FMT", "Transplant"))

4.10.1 Alpha diversity

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("neutral","richness","phylogenetic"))) %>%
  filter(type=="Treatment" & treatment %in% c( "Transplant", "FMT")) %>% 
  ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
  geom_jitter(width = 0.2, show.legend = FALSE) +
  geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
  scale_color_manual(name="treatment",
          breaks=c("FMT","Transplant"),
          values=c('#40714b', "#d5b52c")) +
  scale_fill_manual(name="treatment",
          breaks=c("FMT","Transplant"),
          values=c('#40714b50', "#d5b52c50")) +
  facet_wrap(. ~ metric, scales = "free") +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

Richness

Modelq0GLMMNB <- glmer.nb(richness ~ treatment+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ treatment)

Neutral

Model_neutral <- lme(fixed = neutral ~ treatment, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_neutral)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC     BIC   logLik
  173.0404 178.263 -81.5202

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    4.263247 9.322271

Fixed effects:  neutral ~ treatment 
                         Value Std.Error DF   t-value p-value
(Intercept)          17.310151  3.416951 13  5.065963  0.0002
treatmentAntibiotics -8.868175  4.744329 13 -1.869216  0.0843
treatmentFMT1         7.289358  4.552796 13  1.601073  0.1334
 Correlation: 
                     (Intr) trtmnA
treatmentAntibiotics -0.596       
treatmentFMT1        -0.621  0.457

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.93185401 -0.56384573 -0.04015247  0.47657219  1.55593253 

Number of Observations: 24
Number of Groups: 9 

Phylogenetic

Model_phylo <- lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_phylo)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
      AIC      BIC   logLik
  92.7946 98.01721 -41.3973

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:   0.5507824 1.405346

Fixed effects:  phylogenetic ~ treatment 
                         Value Std.Error DF   t-value p-value
(Intercept)           5.488256 0.5031412 13 10.907984  0.0000
treatmentAntibiotics -1.198294 0.7137109 13 -1.678962  0.1170
treatmentFMT1        -1.351708 0.6855445 13 -1.971729  0.0703
 Correlation: 
                     (Intr) trtmnA
treatmentAntibiotics -0.611       
treatmentFMT1        -0.636  0.456

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.00120146 -0.44426020  0.00785757  0.38719755  1.97397308 

Number of Observations: 24
Number of Groups: 9 

4.10.2 Beta diversity

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.11100 0.111002 11.752    999  0.001 ***
Residuals 28 0.26447 0.009445                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(richness ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(richness))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 1.184578 0.1548451 5.13002 0.001
Residual 28 6.465508 0.8451549 NA NA
Total 29 7.650086 1.0000000 NA NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.04408 0.044082 3.1744    999  0.093 .
Residuals 28 0.38883 0.013887                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(neutral ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(neutral))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 1.147077 0.1608783 5.368223 0.001
Residual 28 5.983012 0.8391217 NA NA
Total 29 7.130089 1.0000000 NA NA

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00007 0.0000689 0.0044    999  0.941
Residuals 28 0.43637 0.0155847                     
adonis2(phylo ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(phylo))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 0.1298187 0.1018776 3.17615 0.001
Residual 28 1.1444432 0.8981224 NA NA
Total 29 1.2742619 1.0000000 NA NA
p0<-richness %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 4) +
  scale_color_manual(name="treatment",
          breaks=c("FMT","Transplant"),
          values=c('#40714b', "#d5b52c")) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )

p1<-neutral %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 4) +
  scale_color_manual(name="treatment",
          breaks=c("FMT","Transplant"),
          values=c('#40714b', "#d5b52c"))+  
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )

p2<-phylo %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 4) +
  scale_color_manual(name="treatment",
          breaks=c("FMT","Transplant"),
          values=c('#40714b', "#d5b52c")) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )
ggarrange(p0, p1, p2, ncol=3, nrow=1, common.legend = TRUE, legend="right")

4.10.3 Structural zeros

FMT_samples <- sample_metadata %>%
  filter(type == "Treatment" & treatment == "FMT") %>% 
  dplyr::select(sample) %>%
  pull()

Transplant_samples <- sample_metadata %>%
  filter(type == "Treatment" & treatment =="Transplant") %>% 
  dplyr::select(sample) %>% pull()

structural_zeros <- genome_counts_filt %>% 
  select(one_of(c("genome",subset_meta$sample))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>% 
   rowwise() %>% #compute for each row (genome)
   mutate(all_zeros_FMT = all(c_across(all_of(FMT_samples)) == 0)) %>% 
   mutate(all_zeros_Transplant = all(c_across(all_of(Transplant_samples)) == 0)) %>% 
   mutate(average_FMT = mean(c_across(all_of(FMT_samples)), na.rm = TRUE)) %>% 
   mutate(average_Transplant = mean(c_across(all_of(Transplant_samples)), na.rm = TRUE)) %>% 
   filter(all_zeros_FMT == TRUE || all_zeros_Transplant==TRUE)  %>% 
   mutate(present = case_when(
      all_zeros_FMT & !all_zeros_Transplant ~ "Transplant",
      !all_zeros_FMT & all_zeros_Transplant ~ "FMT",
      !all_zeros_FMT & !all_zeros_Transplant ~ "None",
      TRUE ~ NA_character_
    )) %>%
   mutate(average = ifelse(present == "FMT", average_FMT, average_Transplant)) %>%
   dplyr::select(genome, present, average) %>%
   left_join(genome_metadata, by=join_by(genome==genome)) %>%
   arrange(present,-average)

Structural zeros in each group

fmt <- structural_zeros %>% 
  filter(present=="FMT") %>% 
  count(phylum, name = "FMT") %>%
  arrange(desc(FMT)) 

fmt_f <- structural_zeros %>% 
  filter(present=="FMT") %>% 
  count(family, name = "FMT") %>%
  arrange(desc(FMT)) 
structural_zeros %>% 
  filter(present=="Transplant") %>% 
  count(phylum, name = "Transplant") %>%
  arrange(desc(Transplant)) %>% 
  full_join(.,fmt, by="phylum" ) %>%
  mutate(across(everything(), ~ replace_na(., 0))) %>%
  as.data.frame() %>% 
  summarise(across(where(is.numeric), ~ sum(.x, na.rm = TRUE)))
  Transplant FMT
1         52  55

Phyla to which the structural zeros belong in each group

structural_zeros %>% 
  filter(present=="Transplant") %>% 
  count(phylum, name = "Transplant") %>%
  arrange(desc(Transplant)) %>% 
  full_join(.,fmt, by="phylum" ) %>%
  mutate(across(everything(), ~ replace_na(., 0))) %>% 
  tt()
phylum Transplant FMT
p__Bacillota_A 20 21
p__Bacillota 16 10
p__Pseudomonadota 6 11
p__Bacteroidota 3 6
p__Desulfobacterota 2 2
p__Bacillota_B 1 0
p__Chlamydiota 1 0
p__Cyanobacteriota 1 0
p__Fusobacteriota 1 0
p__Verrucomicrobiota 1 2
p__Actinomycetota 0 1
p__Bacillota_C 0 1
p__Campylobacterota 0 1

Families to which the structural zeros belong in each group

structural_zeros %>% 
  filter(present=="Transplant") %>% 
  count(family, name = "Transplant") %>%
  arrange(desc(Transplant)) %>% 
  full_join(.,fmt_f, by="family" ) %>%
  mutate(across(everything(), ~ replace_na(., 0))) %>% 
  tt()
family Transplant FMT
f__Lachnospiraceae 7 9
f__Erysipelotrichaceae 6 5
f__UBA660 6 0
f__Enterobacteriaceae 5 2
f__CAG-508 3 0
f__Ruminococcaceae 3 3
f__Anaerovoracaceae 2 0
f__Coprobacillaceae 2 0
f__Desulfovibrionaceae 2 2
f__Oscillospiraceae 2 1
f__Tannerellaceae 2 1
f__UBA1242 2 0
f__ 1 3
f__Akkermansiaceae 1 0
f__CAG-239 1 2
f__Chlamydiaceae 1 0
f__Enterococcaceae 1 3
f__Fusobacteriaceae 1 0
f__Gastranaerophilaceae 1 0
f__Marinifilaceae 1 0
f__Mycoplasmoidaceae 1 1
f__Peptococcaceae 1 0
f__Anaerotignaceae 0 2
f__Bacteroidaceae 0 2
f__LL51 0 2
f__UBA3700 0 2
f__Acutalibacteraceae 0 1
f__Arcobacteraceae 0 1
f__CAG-274 0 1
f__CALVMC01 0 1
f__Devosiaceae 0 1
f__Mycobacteriaceae 0 1
f__Pumilibacteraceae 0 1
f__RUG11792 0 1
f__Rhizobiaceae 0 1
f__Rikenellaceae 0 1
f__Sphingobacteriaceae 0 1
f__Streptococcaceae 0 1
f__UBA1997 0 1
f__UBA3830 0 1
f__Weeksellaceae 0 1

4.10.3.1 Functional capacities of the structural zeros

#Aggregate bundle-level GIFTs into the compound level
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% structural_zeros$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>% 
  select_if(~ !is.numeric(.) || sum(.) != 0)

#Aggregate element-level GIFTs into the function level
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)

#Aggregate function-level GIFTs into overall Biosynthesis, Degradation and Structural GIFTs
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)

sample_sub <- sample_metadata %>%
    mutate(treatment=case_when(
    treatment %in% c("Post-FMT1", "Post-FMT2") ~ "FMT",
    TRUE ~ treatment
  ))%>%
  filter(type == "Treatment" & treatment %in% c("FMT", "Transplant"))

genome_counts_row <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% 
  filter(genome %in% structural_zeros$genome) %>% 
  select(one_of(c("genome",sample_sub$sample))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>% 
  column_to_rownames(., "genome") 
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)
element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))
uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

significant_elements <- element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)%>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "sample")%>% 
  left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))

difference_table <- element_gift_filt %>%
  select(-sample) %>%
  group_by(treatment) %>%
  summarise(across(everything(), mean)) %>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric) %>%
  rownames_to_column(., "Elements") %>%
  left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>% 
  arrange(Function) %>% 
  mutate(Difference=FMT-Transplant)%>% 
  mutate(group_color = ifelse(Difference <0, "Transplant","FMT")) 
difference_table %>%
  ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) + 
  geom_col() +
#  geom_point(size=4) + 
  #scale_fill_manual() +
  geom_hline(yintercept=0) + 
  coord_flip()+
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        legend.position = "right", 
        panel.background = element_blank(),
          panel.grid.major = element_line(size = 0.15, linetype = 'solid',
                                colour = "grey"))+
  xlab("Microbial Functional Capacity") + 
  ylab("Mean difference")+
  labs(fill="Elevation")

4.10.4 Differential abundance analysis: composition

phylo_samples <- sample_metadata %>%
    mutate(treatment=case_when(
    treatment %in% c("FMT1", "FMT2") ~ "FMT",
    TRUE ~ treatment
  ))%>%
  filter(type == "Treatment" & treatment %in% c("FMT", "Transplant")) %>% 
  column_to_rownames("sample") %>%
  sample_data()
phylo_genome <- genome_counts_filt %>% 
  filter(!genome %in% structural_zeros$genome) %>% 
  select(one_of(c("genome",rownames(phylo_samples)))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
  column_to_rownames("genome") %>% 
  otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>% 
  filter(genome %in% rownames(phylo_genome)) %>% 
  column_to_rownames("genome") %>% 
  dplyr::select(domain,phylum,class,order,family,genus,species) %>% 
  as.matrix() %>% 
  tax_table()

physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)
ancom_rand_output = ancombc2(data = physeq_genome_filtered, 
                  assay_name = "counts", 
                  tax_level = NULL,
                  fix_formula = "treatment",
                  p_adj_method = "holm", 
                  pseudo_sens = TRUE,
                  prv_cut =0, 
                  lib_cut = 0, 
                  s0_perc = 0.05,
                  group = NULL, 
                  struc_zero = FALSE, 
                  neg_lb = FALSE,
                  alpha = 0.05, 
                  n_cl = 2, 
                  verbose = TRUE,
                  global = FALSE, 
                  pairwise = FALSE, 
                  dunnet = FALSE, 
                  trend = FALSE,
                  iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
                  em_control = list(tol = 1e-5, max_iter = 100),
                  lme_control = lme4::lmerControl(),
                  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100), 
                  trend_control = NULL)
save(ancom_rand_output, file = "data/ancombc_FMT_transplant.Rdata")
load("data/ancombc_FMT_transplant.Rdata")
taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
  rownames_to_column(., "taxon") %>%
  mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))

ancombc_rand_table_mag <- ancom_rand_output$res %>%
  dplyr::select(taxon, lfc_treatmentTransplant, p_treatmentTransplant) %>%
  filter(p_treatmentTransplant < 0.05) %>%
  dplyr::arrange(p_treatmentTransplant) %>%
  merge(., taxonomy, by="taxon") %>%
  mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
  dplyr::arrange(lfc_treatmentTransplant)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", ""))  %>%
  right_join(taxonomy, by=join_by(phylum == phylum)) %>%
  dplyr::select(phylum, colors) %>%
  mutate(colors = str_c(colors, "80"))  %>% #add 80% alpha
    unique() %>%
    dplyr::arrange(phylum)

tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
  
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
    dplyr::arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

Differential abundance MAGs in each group

fmt_count <- ancombc_rand_table_mag %>%
  filter(lfc_treatmentTransplant<0)  %>% 
  count(phylum, name = "FMT") %>%
  arrange(desc(FMT)) 

ancombc_rand_table_mag %>%
  filter(lfc_treatmentTransplant>0)  %>% 
  count(phylum, name = "Transplant") %>%
  arrange(desc(Transplant))  %>% 
  full_join(.,fmt_count, by="phylum")%>%
  mutate(across(where(is.numeric), ~ replace_na(., 0)))
             phylum Transplant FMT
1      Bacteroidota         13   5
2       Bacillota_A          4  17
3         Bacillota          3   1
4  Campylobacterota          1   1
5  Desulfobacterota          1   2
6     Spirochaetota          1   0
7    Pseudomonadota          0   5
8   Cyanobacteriota          0   2
9       Bacillota_B          0   1
10      Bacillota_C          0   1
ancombc_rand_table_mag %>%
  filter(lfc_treatmentTransplant>0)  %>% 
  count(phylum, name = "Transplant") %>%
  arrange(desc(Transplant))  %>% 
  full_join(.,fmt_count, by="phylum")%>%
  mutate(across(where(is.numeric), ~ replace_na(., 0))) %>% 
  as.data.frame() %>% 
  summarise(across(where(is.numeric), ~ sum(.x, na.rm = TRUE)))
  Transplant FMT
1         23  35
ancombc_rand_table_mag%>%
      mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_treatmentTransplant, y=forcats::fct_reorder(genome,lfc_treatmentTransplant), fill=phylum)) + #forcats::fct_rev()
  geom_col() + 
  scale_fill_manual(values=tax_color) + 
  geom_hline(yintercept=0) + 
#  coord_flip()+
  theme(panel.background = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 8),
        axis.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14, face = "bold"),
        legend.position = "right", legend.box = "vertical")+
  xlab("log2FoldChange") + 
  ylab("Species")+
  guides(fill=guide_legend(title="Phylum"))

4.10.5 Differences in the functional capacities

GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>% 
  select_if(~ !is.numeric(.) || sum(.) != 0)
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)

sample_sub <- sample_metadata %>%
    mutate(treatment=case_when(
    treatment %in% c("Post-FMT1", "Post-FMT2") ~ "FMT",
    TRUE ~ treatment
  ))%>%
  filter(type == "Treatment" & treatment %in% c("FMT", "Transplant"))

genome_counts_row <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% 
    select(one_of(c("genome",sample_sub$sample))) %>%
  column_to_rownames(., "genome") 

GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)

4.10.5.1 Community elements differences:

MCI

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))
# A tibble: 2 × 3
  treatment    MCI     sd
  <chr>      <dbl>  <dbl>
1 FMT        0.359 0.0259
2 Transplant 0.356 0.0432
MCI <- GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.94871, p-value = 0.1561
wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 128, p-value = 0.4826
alternative hypothesis: true location shift is not equal to 0
element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))
uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

significant_elements <- element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)%>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "sample")%>% 
  left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))

difference_table <- element_gift_filt %>%
  select(-sample) %>%
  group_by(treatment) %>%
  summarise(across(everything(), mean)) %>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric) %>%
  rownames_to_column(., "Elements") %>%
  left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>% 
  arrange(Function) %>% 
  mutate(Difference=FMT-Transplant)%>% 
  mutate(group_color = ifelse(Difference <0, "Transplant","FMT")) 
means_gift <- element_gift_filt %>% 
  select(-treatment) %>% 
  pivot_longer(!sample, names_to = "Elements", values_to = "abundance") %>% 
  left_join(sample_metadata, by=join_by(sample==sample)) %>% 
  group_by(treatment, Elements) %>%
  summarise(mean=mean(abundance))

log_fold <- means_gift %>%
  group_by(Elements) %>%
  summarise(
    logfc_fmt_transplant = log2(mean[treatment == "FMT"] / mean[treatment == "Transplant"])
    ) %>% 
  left_join(., difference_table, by="Elements")
difference_table %>%
  ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) + 
  geom_col() +
#  geom_point(size=4) + 
  scale_fill_manual(values=c('#40714b', "#d5b52c")) +
  geom_hline(yintercept=0) + 
  coord_flip()+
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        legend.position = "right", 
        panel.background = element_blank(),
          panel.grid.major = element_line(size = 0.15, linetype = 'solid',
                                colour = "grey"))+
  xlab("Microbial Functional Capacity") + 
  ylab("Mean difference")+
  labs(fill="Treatment")

log_fold %>%
  filter(Elements!="D0611") %>% 
  ggplot(aes(x=forcats::fct_reorder(Function,logfc_fmt_transplant), y=logfc_fmt_transplant, fill=group_color)) + 
  geom_col() +
  #scale_fill_manual() +
  geom_hline(yintercept=0) + 
  coord_flip()+
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        legend.position = "right", 
        panel.background = element_blank(),
          panel.grid.major = element_line(size = 0.15, linetype = 'solid',
                                colour = "grey"))+
  xlab("Microbial Functional Capacity") + 
  ylab("Log_fold")+
  labs(fill="Treatment")

#### Community functions differences

MCI

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))
# A tibble: 2 × 3
  treatment    MCI     sd
  <chr>      <dbl>  <dbl>
1 FMT        0.349 0.0221
2 Transplant 0.355 0.0377
MCI <- GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.87044, p-value = 0.001713
wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 121, p-value = 0.6801
alternative hypothesis: true location shift is not equal to 0
function_gift <- GIFTs_functions_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  merge(., sample_metadata[c(12,13)], by="sample")
unique_funct_db<- GIFT_db[c(3,4,5)] %>% 
  distinct(Code_function, .keep_all = TRUE)

significant_functional <- function_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)%>%
  left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))

significant_functional
# A tibble: 3 × 4
  trait  p_value p_adjust Function               
  <chr>    <dbl>    <dbl> <chr>                  
1 B04   0.000220  0.00441 SCFA biosynthesis      
2 B10   0.00284   0.0189  Antibiotic biosynthesis
3 S02   0.00174   0.0174  Appendages